Merge remote-tracking branch 'origin/master' into complex-strumpack-dev

# Conflicts:
#	examples/ex11p.cpp
#	linalg/strumpack.cpp
#	linalg/strumpack.hpp

examples/CMakeLists.txt

@@ -27,6 +27,7 @@ list(APPEND ALL_EXE_SRCS
   ex18.cpp
   ex19.cpp
   ex20.cpp
+  ex21.cpp
   ex22.cpp
   )
 
@@ -52,6 +53,7 @@ if (MFEM_USE_MPI)
     ex18p.cpp
     ex19p.cpp
     ex20p.cpp
+    ex21p.cpp
     ex22p.cpp
     )
 endif()

examples/ex21.cpp

@@ -0,0 +1,477 @@
+//                               MFEM Example 21
+//
+// Compile with: make ex21
+//
+// Sample runs:  ex21 -m ../data/inline-segment.mesh -o 3
+//               ex21 -m ../data/inline-tri.mesh -o 3
+//               ex21 -m ../data/inline-quad.mesh -o 3
+//               ex21 -m ../data/inline-quad.mesh -o 3 -p 1
+//               ex21 -m ../data/inline-quad.mesh -o 3 -p 2
+//               ex21 -m ../data/inline-tet.mesh -o 2
+//               ex21 -m ../data/inline-hex.mesh -o 2
+//               ex21 -m ../data/inline-hex.mesh -o 2 -p 1
+//               ex21 -m ../data/inline-hex.mesh -o 2 -p 2
+//               ex21 -m ../data/star.mesh -o 2 -sigma 10.0
+//
+// Description:  This example code demonstrates the use of MFEM to define and
+//               solve simple complex-valued linear systems.  We implement three
+//               variants of a damped harmonic oscillator:
+//
+//               1) A scalar H1 field
+//                  -Div(a Grad u) - omega^2 b u + i omega c u = 0
+//
+//               2) A vector H(Curl) field
+//                  Curl(a Curl u) - omega^2 b u + i omega c u = 0
+//
+//               3) A vector H(Div) field
+//                  -Grad(a Div u) - omega^2 b u + i omega c u = 0
+//
+//               In each case the field is driven by a forced oscillation, with
+//               angular frequency omega, imposed at the boundary or a portion
+//               of the boundary.
+//
+//               In electromagnetics the coefficients are typically named the
+//               permeability, mu = 1/a, permittivity, epsilon = b, and
+//               conductivity, sigma = c.  The user can specify these constants
+//               using either set of names.
+//
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+static double mu_ = 1.0;
+static double epsilon_ = 1.0;
+static double sigma_ = 20.0;
+static double omega_ = 10.0;
+
+double u0_real_exact(const Vector &);
+double u0_imag_exact(const Vector &);
+
+void u1_real_exact(const Vector &, Vector &);
+void u1_imag_exact(const Vector &, Vector &);
+
+void u2_real_exact(const Vector &, Vector &);
+void u2_imag_exact(const Vector &, Vector &);
+
+bool check_for_inline_mesh(const char * mesh_file);
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   const char *mesh_file = "../data/inline-quad.mesh";
+   int ref_levels = 0;
+   int order = 1;
+   int prob = 0;
+   double freq = -1.0;
+   double a_coef = 0.0;
+   bool visualization = 1;
+   bool herm_conv = true;
+   bool exact_sol = true;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&ref_levels, "-r", "--refine",
+                  "Number of times to refine the mesh uniformly.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree).");
+   args.AddOption(&prob, "-p", "--problem-type",
+                  "Choose from 0: H_1, 1: H(Curl), or 2: H(Div) "
+                  "damped harmonic oscillator.");
+   args.AddOption(&a_coef, "-a", "--stiffness-coef",
+                  "Stiffness coefficient (spring constant or 1/mu).");
+   args.AddOption(&epsilon_, "-b", "--mass-coef",
+                  "Mass coefficient (or epsilon).");
+   args.AddOption(&sigma_, "-c", "--damping-coef",
+                  "Damping coefficient (or sigma).");
+   args.AddOption(&mu_, "-mu", "--permeability",
+                  "Permeability of free space (or 1/(spring constant)).");
+   args.AddOption(&epsilon_, "-eps", "--permittivity",
+                  "Permittivity of free space (or mass constant).");
+   args.AddOption(&sigma_, "-sigma", "--conductivity",
+                  "Conductivity (or damping constant).");
+   args.AddOption(&freq, "-f", "--frequency",
+                  "Frequency (in Hz).");
+   args.AddOption(&herm_conv, "-herm", "--hermitian", "-no-herm",
+                  "--no-hermitian", "Use convention for Hermitian operators.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   if ( a_coef != 0.0 )
+   {
+      mu_ = 1.0 / a_coef;
+   }
+   if ( freq > 0.0 )
+   {
+      omega_ = 2.0 * M_PI * freq;
+   }
+
+   exact_sol = check_for_inline_mesh(mesh_file);
+   if (exact_sol)
+   {
+      cout << "Identified an 'inline' mesh" << endl;
+   }
+
+   ComplexOperator::Convention conv =
+      herm_conv ? ComplexOperator::HERMITIAN : ComplexOperator::BLOCK_SYMMETRIC;
+
+   // 2. Read the mesh from the given mesh file. We can handle triangular,
+   //    quadrilateral, tetrahedral, hexahedral, surface and volume meshes
+   //    with the same code.
+   Mesh *mesh = new Mesh(mesh_file, 1, 1);
+   int dim = mesh->Dimension();
+
+   // 3. Refine the mesh to increase resolution. In this example we do
+   //    'ref_levels' of uniform refinement where the user specifies
+   //    the number of levels with the '-r' option.
+   for (int l = 0; l < ref_levels; l++)
+   {
+      mesh->UniformRefinement();
+   }
+
+   // 4. Define a finite element space on the mesh. Here we use continuous
+   //    Lagrange, Nedelec, or Raviart-Thomas finite elements of the specified
+   //    order.
+   if (dim == 1 && prob != 0 )
+   {
+      cout << "Switching to problem type 0, H1 basis functions, "
+           << "for 1 dimensional mesh." << endl;
+      prob = 0;
+   }
+
+   FiniteElementCollection *fec;
+   switch (prob)
+   {
+      case 0:  fec = new H1_FECollection(order, dim);      break;
+      case 1:  fec = new ND_FECollection(order, dim);      break;
+      case 2:  fec = new RT_FECollection(order - 1, dim);  break;
+   }
+   FiniteElementSpace *fespace = new FiniteElementSpace(mesh, fec);
+   cout << "Number of finite element unknowns: " << fespace->GetTrueVSize()
+        << endl;
+
+   // 5. Determine the list of true (i.e. conforming) essential boundary dofs.
+   //    In this example, the boundary conditions are defined based on the type
+   //    of mesh and the problem type.
+   Array<int> ess_tdof_list;
+   Array<int> ess_bdr;
+   if (mesh->bdr_attributes.Size())
+   {
+      ess_bdr.SetSize(mesh->bdr_attributes.Max());
+      ess_bdr = 1;
+      if (exact_sol)
+      {
+         switch (prob)
+         {
+            case 0:   ess_bdr = 0; ess_bdr[0] = 1;  break;
+            default:  ess_bdr = 1; ess_bdr[2] = 0;  break;
+         }
+      }
+      fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+
+   // 6. Set up the linear form b(.) which corresponds to the
+   //    right-hand side of the FEM linear system.
+   ComplexLinearForm b(fespace, conv);
+   b.Vector::operator=(0.0);
+
+   // 7. Define the solution vector u as a finite element grid function
+   //    corresponding to fespace. Initialize u with initial guess of 1+0i
+   //    or the exact solution if it is known.
+   ComplexGridFunction u(fespace);
+   ComplexGridFunction * u_exact = NULL;
+   if (exact_sol) { u_exact = new ComplexGridFunction(fespace); }
+
+   FunctionCoefficient u0_r(u0_real_exact);
+   FunctionCoefficient u0_i(u0_imag_exact);
+   VectorFunctionCoefficient u1_r(dim, u1_real_exact);
+   VectorFunctionCoefficient u1_i(dim, u1_imag_exact);
+   VectorFunctionCoefficient u2_r(dim, u2_real_exact);
+   VectorFunctionCoefficient u2_i(dim, u2_imag_exact);
+
+   ConstantCoefficient zeroCoef(0.0);
+   ConstantCoefficient oneCoef(1.0);
+
+   Vector zeroVec(dim); zeroVec = 0.0;
+   Vector  oneVec(dim);  oneVec = 0.0; oneVec[(prob==2)?(dim-1):0] = 1.0;
+   VectorConstantCoefficient zeroVecCoef(zeroVec);
+   VectorConstantCoefficient oneVecCoef(oneVec);
+
+   switch (prob)
+   {
+      case 0:
+         u.ProjectBdrCoefficient(oneCoef, zeroCoef, ess_bdr);
+         if (exact_sol) { u_exact->ProjectCoefficient(u0_r, u0_i); }
+         break;
+      case 1:
+         u.ProjectBdrCoefficientTangent(oneVecCoef, zeroVecCoef, ess_bdr);
+         if (exact_sol) { u_exact->ProjectCoefficient(u1_r, u1_i); }
+         break;
+      case 2:
+         u.ProjectBdrCoefficientNormal(oneVecCoef, zeroVecCoef, ess_bdr);
+         if (exact_sol) { u_exact->ProjectCoefficient(u2_r, u2_i); }
+         break;
+   }
+
+   if (visualization && exact_sol)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream sol_sock_r(vishost, visport);
+      socketstream sol_sock_i(vishost, visport);
+      sol_sock_r.precision(8);
+      sol_sock_i.precision(8);
+      sol_sock_r << "solution\n" << *mesh << u_exact->real()
+                 << "window_title 'Exact Real Part'" << flush;
+      sol_sock_i << "solution\n" << *mesh << u_exact->imag()
+                 << "window_title 'Exact Imaginary Part'" << flush;
+   }
+
+   // 8. Set up the sesquilinear form a(.,.) on the finite element
+   //    space corresponding to the damped harmonic oscillator operator
+   //    of the appropriate type:
+   //
+   //    0) A scalar H1 field
+   //       -Div(a Grad) - omega^2 b + i omega c
+   //
+   //    1) A vector H(Curl) field
+   //       Curl(a Curl) - omega^2 b + i omega c
+   //
+   //    2) A vector H(Div) field
+   //       -Grad(a Div) - omega^2 b + i omega c
+   //
+   ConstantCoefficient stiffnessCoef(1.0/mu_);
+   ConstantCoefficient massCoef(-omega_ * omega_ * epsilon_);
+   ConstantCoefficient lossCoef(omega_ * sigma_);
+   ConstantCoefficient negMassCoef(omega_ * omega_ * epsilon_);
+
+   SesquilinearForm *a = new SesquilinearForm(fespace, conv);
+   switch (prob)
+   {
+      case 0:
+         a->AddDomainIntegrator(new DiffusionIntegrator(stiffnessCoef),
+                                NULL);
+         a->AddDomainIntegrator(new MassIntegrator(massCoef),
+                                new MassIntegrator(lossCoef));
+         break;
+      case 1:
+         a->AddDomainIntegrator(new CurlCurlIntegrator(stiffnessCoef),
+                                NULL);
+         a->AddDomainIntegrator(new VectorFEMassIntegrator(massCoef),
+                                new VectorFEMassIntegrator(lossCoef));
+         break;
+      case 2:
+         a->AddDomainIntegrator(new DivDivIntegrator(stiffnessCoef),
+                                NULL);
+         a->AddDomainIntegrator(new VectorFEMassIntegrator(massCoef),
+                                new VectorFEMassIntegrator(lossCoef));
+         break;
+   }
+
+   // 9. Assemble the bilinear form and the corresponding linear
+   //    system, applying any necessary transformations such as:
+   //    assembly, eliminating boundary conditions, applying conforming
+   //    constraints for non-conforming AMR, etc.
+   a->Assemble();
+
+   OperatorHandle A;
+   Vector B, U;
+
+   a->FormLinearSystem(ess_tdof_list, u, b, A, U, B);
+   u = 0.0;
+   U = 0.0;
+
+   {
+      ComplexSparseMatrix * Asp =
+         dynamic_cast<ComplexSparseMatrix*>(A.Ptr());
+
+      cout << "Size of linear system: "
+           << 2 * Asp->real().Width() << endl << endl;
+   }
+
+   // 10. Define and apply a GMRES solver for AU=B.
+   {
+      GMRESSolver gmres;
+      gmres.SetOperator(*A.Ptr());
+      gmres.SetRelTol(1e-12);
+      gmres.SetMaxIter(1000);
+      gmres.SetPrintLevel(1);
+      gmres.Mult(B, U);
+   }
+
+   // 11. Recover the solution as a finite element grid function and
+   //     compute the errors if the exact solution is known.
+   a->RecoverFEMSolution(U, b, u);
+
+   if (exact_sol)
+   {
+      double err_r = -1.0;
+      double err_i = -1.0;
+
+      switch (prob)
+      {
+         case 0:
+            err_r = u.real().ComputeL2Error(u0_r);
+            err_i = u.imag().ComputeL2Error(u0_i);
+            break;
+         case 1:
+            err_r = u.real().ComputeL2Error(u1_r);
+            err_i = u.imag().ComputeL2Error(u1_i);
+            break;
+         case 2:
+            err_r = u.real().ComputeL2Error(u2_r);
+            err_i = u.imag().ComputeL2Error(u2_i);
+            break;
+      }
+
+      cout << endl;
+      cout << "|| Re (u_h - u) ||_{L^2} = " << err_r << endl;
+      cout << "|| Im (u_h - u) ||_{L^2} = " << err_i << endl;
+      cout << endl;
+   }
+
+   // 12. Save the refined mesh and the solution. This output can be
+   //     viewed later using GLVis: "glvis -m mesh -g sol".
+   {
+      ofstream mesh_ofs("refined.mesh");
+      mesh_ofs.precision(8);
+      mesh->Print(mesh_ofs);
+
+      ofstream sol_r_ofs("sol_r.gf");
+      ofstream sol_i_ofs("sol_i.gf");
+      sol_r_ofs.precision(8);
+      sol_i_ofs.precision(8);
+      u.real().Save(sol_r_ofs);
+      u.imag().Save(sol_i_ofs);
+   }
+
+   // 13. Send the solution by socket to a GLVis server.
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream sol_sock_r(vishost, visport);
+      socketstream sol_sock_i(vishost, visport);
+      sol_sock_r.precision(8);
+      sol_sock_i.precision(8);
+      sol_sock_r << "solution\n" << *mesh << u.real()
+                 << "window_title 'Comp Real Part'" << flush;
+      sol_sock_i << "solution\n" << *mesh << u.imag()
+                 << "window_title 'Comp Imaginary Part'" << flush;
+   }
+   if (visualization && exact_sol)
+   {
+      *u_exact -= u;
+
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream sol_sock_r(vishost, visport);
+      socketstream sol_sock_i(vishost, visport);
+      sol_sock_r.precision(8);
+      sol_sock_i.precision(8);
+      sol_sock_r << "solution\n" << *mesh << u_exact->real()
+                 << "window_title 'Exact-Comp Real Part'" << flush;
+      sol_sock_i << "solution\n" << *mesh << u_exact->imag()
+                 << "window_title 'Exact-Comp Imaginary Part'" << flush;
+   }
+   if (visualization)
+   {
+      GridFunction u_t(fespace);
+      u_t = u.real();
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream sol_sock(vishost, visport);
+      sol_sock.precision(8);
+      sol_sock << "solution\n" << *mesh << u_t
+               << "window_title 'Harmonic Solution (t = 0.0 T)'"
+               << "pause\n" << flush;
+
+      cout << "GLVis visualization paused."
+           << " Press space (in the GLVis window) to resume it.\n";
+      int num_frames = 32;
+      int i = 0;
+      while (sol_sock)
+      {
+         double t = (double)(i % num_frames) / num_frames;
+         ostringstream oss;
+         oss << "Harmonic Solution (t = " << t << " T)";
+
+         add(cos( 2.0 * M_PI * t), u.real(),
+             sin(-2.0 * M_PI * t), u.imag(), u_t);
+         sol_sock << "solution\n" << *mesh << u_t
+                  << "window_title '" << oss.str() << "'" << flush;
+         i++;
+      }
+   }
+
+   // 14. Free the used memory.
+   delete a;
+   delete u_exact;
+   delete fespace;
+   delete fec;
+   delete mesh;
+
+   return 0;
+}
+
+bool check_for_inline_mesh(const char * mesh_file)
+{
+   string file(mesh_file);
+   size_t p0 = file.find_last_of("/");
+   string s0 = file.substr((p0==string::npos)?0:(p0+1),7);
+   return s0 == "inline-";
+}
+
+complex<double> u0_exact(const Vector &x)
+{
+   int dim = x.Size();
+   complex<double> i(0.0, 1.0);
+   complex<double> alpha = (epsilon_ * omega_ - i * sigma_);
+   complex<double> kappa = std::sqrt(mu_ * omega_* alpha);
+   return std::exp(-i * kappa * x[dim - 1]);
+}
+
+double u0_real_exact(const Vector &x)
+{
+   return u0_exact(x).real();
+}
+
+double u0_imag_exact(const Vector &x)
+{
+   return u0_exact(x).imag();
+}
+
+void u1_real_exact(const Vector &x, Vector &v)
+{
+   int dim = x.Size();
+   v.SetSize(dim); v = 0.0; v[0] = u0_real_exact(x);
+}
+
+void u1_imag_exact(const Vector &x, Vector &v)
+{
+   int dim = x.Size();
+   v.SetSize(dim); v = 0.0; v[0] = u0_imag_exact(x);
+}
+
+void u2_real_exact(const Vector &x, Vector &v)
+{
+   int dim = x.Size();
+   v.SetSize(dim); v = 0.0; v[dim-1] = u0_real_exact(x);
+}
+
+void u2_imag_exact(const Vector &x, Vector &v)
+{
+   int dim = x.Size();
+   v.SetSize(dim); v = 0.0; v[dim-1] = u0_imag_exact(x);
+}

examples/ex21p.cpp

@@ -0,0 +1,658 @@
+//                       MFEM Example 21 - Parallel Version
+//
+// Compile with: make ex21p
+//
+// Sample runs:  mpirun -np 4 ex21p -m ../data/inline-segment.mesh -o 3
+//               mpirun -np 4 ex21p -m ../data/inline-tri.mesh -o 3
+//               mpirun -np 4 ex21p -m ../data/inline-quad.mesh -o 3
+//               mpirun -np 4 ex21p -m ../data/inline-quad.mesh -o 3 -p 1
+//               mpirun -np 4 ex21p -m ../data/inline-quad.mesh -o 3 -p 2
+//               mpirun -np 4 ex21p -m ../data/inline-tet.mesh -o 2
+//               mpirun -np 4 ex21p -m ../data/inline-hex.mesh -o 2
+//               mpirun -np 4 ex21p -m ../data/inline-hex.mesh -o 2 -p 1
+//               mpirun -np 4 ex21p -m ../data/inline-hex.mesh -o 2 -p 2
+//               mpirun -np 4 ex21p -m ../data/star.mesh -o 2 -sigma 10.0
+//
+// Description:  This example code demonstrates the use of MFEM to define and
+//               solve simple complex-valued linear systems.  We implement three
+//               variants of a damped harmonic oscillator:
+//
+//               1) A scalar H1 field
+//                  -Div(a Grad u) - omega^2 b u + i omega c u = 0
+//
+//               2) A vector H(Curl) field
+//                  Curl(a Curl u) - omega^2 b u + i omega c u = 0
+//
+//               3) A vector H(Div) field
+//                  -Grad(a Div u) - omega^2 b u + i omega c u = 0
+//
+//               In each case the field is driven by a forced oscillation, with
+//               angular frequency omega, imposed at the boundary or a portion
+//               of the boundary.
+//
+//               In electromagnetics the coefficients are typically named the
+//               permeability, mu = 1/a, permittivity, epsilon = b, and
+//               conductivity, sigma = c.  The user can specify these constants
+//               using either set of names.
+//
+//#define MFEM_STRUMPACK_SRC
+
+#include <fstream>
+#include <iostream>
+#include "mfem.hpp"
+
+using namespace std;
+using namespace mfem;
+
+static double mu_ = 1.0;
+static double epsilon_ = 1.0;
+static double sigma_ = 20.0;
+static double omega_ = 10.0;
+
+double u0_real_exact(const Vector &);
+double u0_imag_exact(const Vector &);
+
+void u1_real_exact(const Vector &, Vector &);
+void u1_imag_exact(const Vector &, Vector &);
+
+void u2_real_exact(const Vector &, Vector &);
+void u2_imag_exact(const Vector &, Vector &);
+
+bool check_for_inline_mesh(const char * mesh_file);
+
+int main(int argc, char *argv[])
+{
+   // 1. Initialize MPI.
+   int num_procs, myid;
+   MPI_Init(&argc, &argv);
+   MPI_Comm comm = MPI_COMM_WORLD;
+   MPI_Comm_size(comm, &num_procs);
+   MPI_Comm_rank(comm, &myid);
+
+   // 2. Parse command-line options.
+   const char *mesh_file = "../data/inline-quad.mesh";
+   int ser_ref_levels = 1;
+   int par_ref_levels = 1;
+   int order = 1;
+   int prob = 0;
+   double freq = -1.0;
+   double a_coef = 0.0;
+   bool visualization = 1;
+   bool herm_conv = true;
+   bool exact_sol = true;
+#ifdef MFEM_USE_STRUMPACK
+   bool strumpack = false;
+#endif
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
+                  "Number of times to refine the mesh uniformly in serial.");
+   args.AddOption(&par_ref_levels, "-rp", "--refine-parallel",
+                  "Number of times to refine the mesh uniformly in parallel.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree).");
+   args.AddOption(&prob, "-p", "--problem-type",
+                  "Choose from 0: H_1, 1: H(Curl), or 2: H(Div) "
+                  "damped harmonic oscillator.");
+   args.AddOption(&a_coef, "-a", "--stiffness-coef",
+                  "Stiffness coefficient (spring constant or 1/mu).");
+   args.AddOption(&epsilon_, "-b", "--mass-coef",
+                  "Mass coefficient (or epsilon).");
+   args.AddOption(&sigma_, "-c", "--damping-coef",
+                  "Damping coefficient (or sigma).");
+   args.AddOption(&mu_, "-mu", "--permeability",
+                  "Permeability of free space (or 1/(spring constant)).");
+   args.AddOption(&epsilon_, "-eps", "--permittivity",
+                  "Permittivity of free space (or mass constant).");
+   args.AddOption(&sigma_, "-sigma", "--conductivity",
+                  "Conductivity (or damping constant).");
+   args.AddOption(&freq, "-f", "--frequency",
+                  "Frequency (in Hz).");
+#ifdef MFEM_USE_STRUMPACK
+   args.AddOption(&strumpack, "-strumpack", "--strumpack-solver",
+                  "-no-strumpack", "--no-strumpack-solver",
+                  "Use STRUMPACK's double complex linear solver.");
+#endif
+   args.AddOption(&herm_conv, "-herm", "--hermitian", "-no-herm",
+                  "--no-hermitian", "Use convention for Hermitian operators.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      MPI_Finalize();
+      return 1;
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+
+   if ( a_coef != 0.0 )
+   {
+      mu_ = 1.0 / a_coef;
+   }
+   if ( freq > 0.0 )
+   {
+      omega_ = 2.0 * M_PI * freq;
+   }
+
+   exact_sol = check_for_inline_mesh(mesh_file);
+   if (myid == 0 && exact_sol)
+   {
+      cout << "Identified an 'inline' mesh" << endl;
+   }
+
+   ComplexOperator::Convention conv =
+      herm_conv ? ComplexOperator::HERMITIAN : ComplexOperator::BLOCK_SYMMETRIC;
+
+   // 3. Read the (serial) mesh from the given mesh file on all processors.  We
+   //    can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
+   //    and volume meshes with the same code.
+   Mesh *mesh = new Mesh(mesh_file, 1, 1);
+   int dim = mesh->Dimension();
+
+   // 4. Refine the serial mesh on all processors to increase the resolution.
+   for (int l = 0; l < ser_ref_levels; l++)
+   {
+      mesh->UniformRefinement();
+   }
+
+   // 5. Define a parallel mesh by a partitioning of the serial mesh. Refine
+   //    this mesh further in parallel to increase the resolution. Once the
+   //    parallel mesh is defined, the serial mesh can be deleted.
+   ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
+   delete mesh;
+   for (int l = 0; l < par_ref_levels; l++)
+   {
+      pmesh->UniformRefinement();
+   }
+
+   // 6. Define a parallel finite element space on the parallel
+   //    mesh. Here we use continuous Lagrange, Nedelec, or
+   //    Raviart-Thomas finite elements of the specified order.
+   if (dim == 1 && prob != 0 )
+   {
+      if (myid == 0)
+      {
+         cout << "Switching to problem type 0, H1 basis functions, "
+              << "for 1 dimensional mesh." << endl;
+      }
+      prob = 0;
+   }
+
+   FiniteElementCollection *fec;
+   switch (prob)
+   {
+      case 0:  fec = new H1_FECollection(order, dim);      break;
+      case 1:  fec = new ND_FECollection(order, dim);      break;
+      case 2:  fec = new RT_FECollection(order - 1, dim);  break;
+   }
+   ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, fec);
+   HYPRE_Int size = fespace->GlobalTrueVSize();
+   if (myid == 0)
+   {
+      cout << "Number of finite element unknowns: " << size << endl;
+   }
+
+   // 7. Determine the list of true (i.e. parallel conforming) essential
+   //    boundary dofs. In this example, the boundary conditions are defined
+   //    based on the type of mesh and the problem type.
+   Array<int> ess_tdof_list;
+   Array<int> ess_bdr;
+   if (pmesh->bdr_attributes.Size())
+   {
+      ess_bdr.SetSize(pmesh->bdr_attributes.Max());
+      ess_bdr = 1;
+      if (exact_sol)
+      {
+         switch (prob)
+         {
+            case 0:   ess_bdr = 0; ess_bdr[0] = 1;  break;
+            default:  ess_bdr = 1; ess_bdr[2] = 0;  break;
+         }
+      }
+      fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+
+   // 8. Set up the parallel linear form b(.) which corresponds to the
+   //    right-hand side of the FEM linear system.
+   ParComplexLinearForm b(fespace, conv);
+   b.Vector::operator=(0.0);
+
+   // 9. Define the solution vector u as a parallel finite element
+   //    grid function corresponding to fespace. Initialize u with
+   //    initial guess of 1+0i or the exact solution if it is known.
+   ParComplexGridFunction u(fespace);
+   ParComplexGridFunction * u_exact = NULL;
+   if (exact_sol) { u_exact = new ParComplexGridFunction(fespace); }
+
+   FunctionCoefficient u0_r(u0_real_exact);
+   FunctionCoefficient u0_i(u0_imag_exact);
+   VectorFunctionCoefficient u1_r(dim, u1_real_exact);
+   VectorFunctionCoefficient u1_i(dim, u1_imag_exact);
+   VectorFunctionCoefficient u2_r(dim, u2_real_exact);
+   VectorFunctionCoefficient u2_i(dim, u2_imag_exact);
+
+   ConstantCoefficient zeroCoef(0.0);
+   ConstantCoefficient oneCoef(1.0);
+
+   Vector zeroVec(dim); zeroVec = 0.0;
+   Vector  oneVec(dim);  oneVec = 0.0; oneVec[(prob==2)?(dim-1):0] = 1.0;
+   VectorConstantCoefficient zeroVecCoef(zeroVec);
+   VectorConstantCoefficient oneVecCoef(oneVec);
+
+   switch (prob)
+   {
+      case 0:
+         u.ProjectBdrCoefficient(oneCoef, zeroCoef, ess_bdr);
+         if (exact_sol) { u_exact->ProjectCoefficient(u0_r, u0_i); }
+         break;
+      case 1:
+         u.ProjectBdrCoefficientTangent(oneVecCoef, zeroVecCoef, ess_bdr);
+         if (exact_sol) { u_exact->ProjectCoefficient(u1_r, u1_i); }
+         break;
+      case 2:
+         u.ProjectBdrCoefficientNormal(oneVecCoef, zeroVecCoef, ess_bdr);
+         if (exact_sol) { u_exact->ProjectCoefficient(u2_r, u2_i); }
+         break;
+   }
+
+   if (visualization && exact_sol)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream sol_sock_r(vishost, visport);
+      socketstream sol_sock_i(vishost, visport);
+      sol_sock_r << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock_i << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock_r.precision(8);
+      sol_sock_i.precision(8);
+      sol_sock_r << "solution\n" << *pmesh << u_exact->real()
+                 << "window_title 'Exact Real Part'" << flush;
+      sol_sock_i << "solution\n" << *pmesh << u_exact->imag()
+                 << "window_title 'Exact Imaginary Part'" << flush;
+   }
+
+   // 10. Set up the parallel sesquilinear form a(.,.) on the finite element
+   //     space corresponding to the damped harmonic oscillator operator
+   //     of the appropriate type:
+   //
+   //     0) A scalar H1 field
+   //        -Div(a Grad) - omega^2 b + i omega c
+   //
+   //     1) A vector H(Curl) field
+   //        Curl(a Curl) - omega^2 b + i omega c
+   //
+   //     2) A vector H(Div) field
+   //        -Grad(a Div) - omega^2 b + i omega c
+   //
+   ConstantCoefficient stiffnessCoef(1.0/mu_);
+   ConstantCoefficient massCoef(-omega_ * omega_ * epsilon_);
+   ConstantCoefficient lossCoef(omega_ * sigma_);
+   ConstantCoefficient negMassCoef(omega_ * omega_ * epsilon_);
+
+   ParSesquilinearForm *a = new ParSesquilinearForm(fespace, conv);
+   switch (prob)
+   {
+      case 0:
+         a->AddDomainIntegrator(new DiffusionIntegrator(stiffnessCoef),
+                                NULL);
+         a->AddDomainIntegrator(new MassIntegrator(massCoef),
+                                new MassIntegrator(lossCoef));
+         break;
+      case 1:
+         a->AddDomainIntegrator(new CurlCurlIntegrator(stiffnessCoef),
+                                NULL);
+         a->AddDomainIntegrator(new VectorFEMassIntegrator(massCoef),
+                                new VectorFEMassIntegrator(lossCoef));
+         break;
+      case 2:
+         a->AddDomainIntegrator(new DivDivIntegrator(stiffnessCoef),
+                                NULL);
+         a->AddDomainIntegrator(new VectorFEMassIntegrator(massCoef),
+                                new VectorFEMassIntegrator(lossCoef));
+         break;
+   }
+
+   // 10a. Set up the parallel bilinear form for the preconditioner
+   //      corresponding to the appropriate operator if the STRUMPACK solver
+   //      has not been selected.
+   //
+   //      0) A scalar H1 field
+   //         -Div(a Grad) - omega^2 b + omega c
+   //
+   //      1) A vector H(Curl) field
+   //         Curl(a Curl) + omega^2 b + omega c
+   //
+   //      2) A vector H(Div) field
+   //         -Grad(a Div) - omega^2 b + omega c
+   //
+   ParBilinearForm *pcOp = NULL;
+#ifdef MFEM_USE_STRUMPACK
+   if (!strumpack)
+#endif
+   {
+      pcOp = new ParBilinearForm(fespace);
+      switch (prob)
+      {
+         case 0:
+            pcOp->AddDomainIntegrator(new DiffusionIntegrator(stiffnessCoef));
+            pcOp->AddDomainIntegrator(new MassIntegrator(massCoef));
+            pcOp->AddDomainIntegrator(new MassIntegrator(lossCoef));
+            break;
+         case 1:
+            pcOp->AddDomainIntegrator(new CurlCurlIntegrator(stiffnessCoef));
+            pcOp->AddDomainIntegrator(new VectorFEMassIntegrator(negMassCoef));
+            pcOp->AddDomainIntegrator(new VectorFEMassIntegrator(lossCoef));
+            break;
+         case 2:
+            pcOp->AddDomainIntegrator(new DivDivIntegrator(stiffnessCoef));
+            pcOp->AddDomainIntegrator(new VectorFEMassIntegrator(massCoef));
+            pcOp->AddDomainIntegrator(new VectorFEMassIntegrator(lossCoef));
+            break;
+      }
+   }
+
+   // 11. Assemble the parallel bilinear form and the corresponding linear
+   //     system, applying any necessary transformations such as: parallel
+   //     assembly, eliminating boundary conditions, applying conforming
+   //     constraints for non-conforming AMR, etc.
+   a->Assemble();
+   if (pcOp) { pcOp->Assemble(); }
+
+   OperatorHandle A;
+   Vector B, U;
+
+   a->FormLinearSystem(ess_tdof_list, u, b, A, U, B);
+   u = 0.0;
+   U = 0.0;
+
+   OperatorHandle PCOp;
+   if (pcOp) { pcOp->FormSystemMatrix(ess_tdof_list, PCOp); }
+
+   if (myid == 0)
+   {
+      ComplexHypreParMatrix * Ahyp =
+         dynamic_cast<ComplexHypreParMatrix*>(A.Ptr());
+
+      cout << "Size of linear system: "
+           << 2 * Ahyp->real().GetGlobalNumRows() << endl << endl;
+   }
+
+   // 12. Define and apply a parallel FGMRES solver for AU=B with a
+   //     block diagonal preconditioner based on the appropriate multigrid
+   //     preconditioner from hypre or simply use STRUMPACK.
+#ifdef MFEM_USE_STRUMPACK
+   if (!strumpack)
+#endif
+   {
+      Array<HYPRE_Int> blockTrueOffsets;
+      blockTrueOffsets.SetSize(3);
+      blockTrueOffsets[0] = 0;
+      blockTrueOffsets[1] = PCOp.Ptr()->Height();
+      blockTrueOffsets[2] = PCOp.Ptr()->Height();
+      blockTrueOffsets.PartialSum();
+
+      BlockDiagonalPreconditioner BDP(blockTrueOffsets);
+
+      Operator * pc_r = NULL;
+      Operator * pc_i = NULL;
+
+      switch (prob)
+      {
+         case 0:
+            pc_r =
+               new HypreBoomerAMG(dynamic_cast<HypreParMatrix&>(*PCOp.Ptr()));
+            pc_i = new ScaledOperator(pc_r,
+                                      (conv == ComplexOperator::HERMITIAN) ?
+                                      1.0:-1.0);
+            break;
+         case 1:
+            pc_r = new HypreAMS(dynamic_cast<HypreParMatrix&>(*PCOp.Ptr()),
+                                fespace);
+            pc_i = new ScaledOperator(pc_r,
+                                      (conv == ComplexOperator::HERMITIAN) ?
+                                      1.0:-1.0);
+            break;
+         case 2:
+            if (dim == 2 )
+            {
+               pc_r = new HypreAMS(dynamic_cast<HypreParMatrix&>(*PCOp.Ptr()),
+                                   fespace);
+            }
+            else
+            {
+               pc_r = new HypreADS(dynamic_cast<HypreParMatrix&>(*PCOp.Ptr()),
+                                   fespace);
+            }
+            pc_i = new ScaledOperator(pc_r,
+                                      (conv == ComplexOperator::HERMITIAN) ?
+                                      1.0:-1.0);
+            break;
+      }
+      BDP.SetDiagonalBlock(0, pc_r);
+      BDP.SetDiagonalBlock(1, pc_i);
+      BDP.owns_blocks = 0;
+
+      FGMRESSolver fgmres(MPI_COMM_WORLD);
+      fgmres.SetPreconditioner(BDP);
+      fgmres.SetOperator(*A.Ptr());
+      fgmres.SetRelTol(1e-12);
+      fgmres.SetMaxIter(1000);
+      fgmres.SetPrintLevel(1);
+      fgmres.Mult(B, U);
+   }
+#ifdef MFEM_USE_STRUMPACK
+   else
+   {
+      ComplexHypreParMatrix * Ahyp =
+         dynamic_cast<ComplexHypreParMatrix*>(A.Ptr());
+
+      STRUMPACKRowLocCmplxMatrix A_strmp(Ahyp->real(), Ahyp->imag());
+
+      STRUMPACKCmplxSolver strmp(argc, argv, comm);
+
+      strmp.SetPrintFactorStatistics(true);
+      strmp.SetPrintSolveStatistics(true);
+      // strmp.SetKrylovSolver(strumpack::KrylovSolver::AUTO); // core dump
+      strmp.SetKrylovSolver(strumpack::KrylovSolver::DIRECT); // core dump
+      // strmp.SetKrylovSolver(strumpack::KrylovSolver::REFINE); // core dump
+      // strmp.SetKrylovSolver(strumpack::KrylovSolver::PREC_GMRES); // index out of range asserts from strumpack::DenseMatrix
+      // strmp.SetKrylovSolver(strumpack::KrylovSolver::GMRES); // WORKS
+      strmp.SetReorderingStrategy(strumpack::ReorderingStrategy::METIS);
+      strmp.SetOperator(A_strmp);
+      strmp.SetFromCommandLine();
+      strmp.Mult(B, U);
+   }
+#endif
+
+   // 13. Recover the parallel grid function corresponding to U. This is the
+   //     local finite element solution on each processor.
+   a->RecoverFEMSolution(U, b, u);
+
+   if (exact_sol)
+   {
+      double err_r = -1.0;
+      double err_i = -1.0;
+
+      switch (prob)
+      {
+         case 0:
+            err_r = u.real().ComputeL2Error(u0_r);
+            err_i = u.imag().ComputeL2Error(u0_i);
+            break;
+         case 1:
+            err_r = u.real().ComputeL2Error(u1_r);
+            err_i = u.imag().ComputeL2Error(u1_i);
+            break;
+         case 2:
+            err_r = u.real().ComputeL2Error(u2_r);
+            err_i = u.imag().ComputeL2Error(u2_i);
+            break;
+      }
+
+      if ( myid == 0 )
+      {
+         cout << endl;
+         cout << "|| Re (u_h - u) ||_{L^2} = " << err_r << endl;
+         cout << "|| Im (u_h - u) ||_{L^2} = " << err_i << endl;
+         cout << endl;
+      }
+   }
+
+   // 14. Save the refined mesh and the solution in parallel. This output can be
+   //     viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
+   {
+      ostringstream mesh_name, sol_r_name, sol_i_name;
+      mesh_name << "mesh." << setfill('0') << setw(6) << myid;
+      sol_r_name << "sol_r." << setfill('0') << setw(6) << myid;
+      sol_i_name << "sol_i." << setfill('0') << setw(6) << myid;
+
+      ofstream mesh_ofs(mesh_name.str().c_str());
+      mesh_ofs.precision(8);
+      pmesh->Print(mesh_ofs);
+
+      ofstream sol_r_ofs(sol_r_name.str().c_str());
+      ofstream sol_i_ofs(sol_i_name.str().c_str());
+      sol_r_ofs.precision(8);
+      sol_i_ofs.precision(8);
+      u.real().Save(sol_r_ofs);
+      u.imag().Save(sol_i_ofs);
+   }
+
+   // 15. Send the solution by socket to a GLVis server.
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream sol_sock_r(vishost, visport);
+      socketstream sol_sock_i(vishost, visport);
+      sol_sock_r << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock_i << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock_r.precision(8);
+      sol_sock_i.precision(8);
+      sol_sock_r << "solution\n" << *pmesh << u.real()
+                 << "window_title 'Comp Real Part'" << flush;
+      sol_sock_i << "solution\n" << *pmesh << u.imag()
+                 << "window_title 'Comp Imaginary Part'" << flush;
+   }
+   if (visualization && exact_sol)
+   {
+      *u_exact -= u;
+
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream sol_sock_r(vishost, visport);
+      socketstream sol_sock_i(vishost, visport);
+      sol_sock_r << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock_i << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock_r.precision(8);
+      sol_sock_i.precision(8);
+      sol_sock_r << "solution\n" << *pmesh << u_exact->real()
+                 << "window_title 'Exact-Comp Real Part'" << flush;
+      sol_sock_i << "solution\n" << *pmesh << u_exact->imag()
+                 << "window_title 'Exact-Comp Imaginary Part'" << flush;
+   }
+   if (visualization)
+   {
+      ParGridFunction u_t(fespace);
+      u_t = u.real();
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream sol_sock(vishost, visport);
+      sol_sock << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock.precision(8);
+      sol_sock << "solution\n" << *pmesh << u_t
+               << "window_title 'Harmonic Solution (t = 0.0 T)'"
+               << "pause\n" << flush;
+      if (myid == 0)
+         cout << "GLVis visualization paused."
+              << " Press space (in the GLVis window) to resume it.\n";
+      int num_frames = 32;
+      int i = 0;
+      while (sol_sock)
+      {
+         double t = (double)(i % num_frames) / num_frames;
+         ostringstream oss;
+         oss << "Harmonic Solution (t = " << t << " T)";
+
+         add(cos( 2.0 * M_PI * t), u.real(),
+             sin(-2.0 * M_PI * t), u.imag(), u_t);
+         sol_sock << "parallel " << num_procs << " " << myid << "\n";
+         sol_sock << "solution\n" << *pmesh << u_t
+                  << "window_title '" << oss.str() << "'" << flush;
+         i++;
+      }
+   }
+
+   // 16. Free the used memory.
+   delete a;
+   delete u_exact;
+   delete pcOp;
+   delete fespace;
+   delete fec;
+   delete pmesh;
+
+   MPI_Finalize();
+
+   return 0;
+}
+
+bool check_for_inline_mesh(const char * mesh_file)
+{
+   string file(mesh_file);
+   size_t p0 = file.find_last_of("/");
+   string s0 = file.substr((p0==string::npos)?0:(p0+1),7);
+   return s0 == "inline-";
+}
+
+complex<double> u0_exact(const Vector &x)
+{
+   int dim = x.Size();
+   complex<double> i(0.0, 1.0);
+   complex<double> alpha = (epsilon_ * omega_ - i * sigma_);
+   complex<double> kappa = std::sqrt(mu_ * omega_* alpha);
+   return std::exp(-i * kappa * x[dim - 1]);
+}
+
+double u0_real_exact(const Vector &x)
+{
+   return u0_exact(x).real();
+}
+
+double u0_imag_exact(const Vector &x)
+{
+   return u0_exact(x).imag();
+}
+
+void u1_real_exact(const Vector &x, Vector &v)
+{
+   int dim = x.Size();
+   v.SetSize(dim); v = 0.0; v[0] = u0_real_exact(x);
+}
+
+void u1_imag_exact(const Vector &x, Vector &v)
+{
+   int dim = x.Size();
+   v.SetSize(dim); v = 0.0; v[0] = u0_imag_exact(x);
+}
+
+void u2_real_exact(const Vector &x, Vector &v)
+{
+   int dim = x.Size();
+   v.SetSize(dim); v = 0.0; v[dim-1] = u0_real_exact(x);
+}
+
+void u2_imag_exact(const Vector &x, Vector &v)
+{
+   int dim = x.Size();
+   v.SetSize(dim); v = 0.0; v[dim-1] = u0_imag_exact(x);
+}

examples/ex3p_sp.cpp

@@ -0,0 +1,334 @@
+//                       MFEM Example 3 - Parallel Version
+//
+// Compile with: make ex3p
+//
+// Sample runs:  mpirun -np 4 ex3p -m ../data/star.mesh
+//               mpirun -np 4 ex3p -m ../data/square-disc.mesh -o 2
+//               mpirun -np 4 ex3p -m ../data/beam-tet.mesh
+//               mpirun -np 4 ex3p -m ../data/beam-hex.mesh
+//               mpirun -np 4 ex3p -m ../data/escher.mesh
+//               mpirun -np 4 ex3p -m ../data/escher.mesh -o 2
+//               mpirun -np 4 ex3p -m ../data/fichera.mesh
+//               mpirun -np 4 ex3p -m ../data/fichera-q2.vtk
+//               mpirun -np 4 ex3p -m ../data/fichera-q3.mesh
+//               mpirun -np 4 ex3p -m ../data/square-disc-nurbs.mesh
+//               mpirun -np 4 ex3p -m ../data/beam-hex-nurbs.mesh
+//               mpirun -np 4 ex3p -m ../data/amr-quad.mesh -o 2
+//               mpirun -np 4 ex3p -m ../data/amr-hex.mesh
+//               mpirun -np 4 ex3p -m ../data/star-surf.mesh -o 2
+//               mpirun -np 4 ex3p -m ../data/mobius-strip.mesh -o 2 -f 0.1
+//               mpirun -np 4 ex3p -m ../data/klein-bottle.mesh -o 2 -f 0.1
+//
+// Description:  This example code solves a simple electromagnetic diffusion
+//               problem corresponding to the second order definite Maxwell
+//               equation curl curl E + E = f with boundary condition
+//               E x n = <given tangential field>. Here, we use a given exact
+//               solution E and compute the corresponding r.h.s. f.
+//               We discretize with Nedelec finite elements in 2D or 3D.
+//
+//               The example demonstrates the use of H(curl) finite element
+//               spaces with the curl-curl and the (vector finite element) mass
+//               bilinear form, as well as the computation of discretization
+//               error when the exact solution is known. Static condensation is
+//               also illustrated.
+//
+//               We recommend viewing examples 1-2 before viewing this example.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+// Exact solution, E, and r.h.s., f. See below for implementation.
+void E_exact(const Vector &, Vector &);
+void f_exact(const Vector &, Vector &);
+double freq = 1.0, kappa;
+int dim;
+
+int main(int argc, char *argv[])
+{
+   // 1. Initialize MPI.
+   int num_procs, myid;
+   MPI_Init(&argc, &argv);
+   MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
+   MPI_Comm_rank(MPI_COMM_WORLD, &myid);
+
+   // 2. Parse command-line options.
+   const char *mesh_file = "../data/beam-tet.mesh";
+   int order = 1;
+   bool static_cond = false;
+   bool visualization = 1;
+#ifdef MFEM_USE_STRUMPACK
+   bool use_strumpack = false;
+#endif
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree).");
+   args.AddOption(&freq, "-f", "--frequency", "Set the frequency for the exact"
+                  " solution.");
+   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
+                  "--no-static-condensation", "Enable static condensation.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+#ifdef MFEM_USE_STRUMPACK
+   args.AddOption(&use_strumpack, "-strumpack", "--strumpack-solver",
+                  "-no-strumpack", "--no-strumpack-solver",
+                  "Use STRUMPACK's double complex linear solver.");
+#endif
+
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      MPI_Finalize();
+      return 1;
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+   kappa = freq * M_PI;
+
+   // 3. Read the (serial) mesh from the given mesh file on all processors.  We
+   //    can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
+   //    and volume meshes with the same code.
+   Mesh *mesh = new Mesh(mesh_file, 1, 1);
+   dim = mesh->Dimension();
+   int sdim = mesh->SpaceDimension();
+
+   // 4. Refine the serial mesh on all processors to increase the resolution. In
+   //    this example we do 'ref_levels' of uniform refinement. We choose
+   //    'ref_levels' to be the largest number that gives a final mesh with no
+   //    more than 1,000 elements.
+   {
+      int ref_levels =
+         (int)floor(log(100000./mesh->GetNE())/log(2.)/dim);
+      for (int l = 0; l < ref_levels; l++)
+      {
+         mesh->UniformRefinement();
+      }
+   }
+
+   // 5. Define a parallel mesh by a partitioning of the serial mesh. Refine
+   //    this mesh further in parallel to increase the resolution. Once the
+   //    parallel mesh is defined, the serial mesh can be deleted. Tetrahedral
+   //    meshes need to be reoriented before we can define high-order Nedelec
+   //    spaces on them.
+   ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
+   delete mesh;
+   {
+      int par_ref_levels = 2;
+      for (int l = 0; l < par_ref_levels; l++)
+      {
+         pmesh->UniformRefinement();
+      }
+   }
+   pmesh->ReorientTetMesh();
+
+   // 6. Define a parallel finite element space on the parallel mesh. Here we
+   //    use the Nedelec finite elements of the specified order.
+   FiniteElementCollection *fec = new ND_FECollection(order, dim);
+   ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, fec);
+   HYPRE_Int size = fespace->GlobalTrueVSize();
+   if (myid == 0)
+   {
+      cout << "Number of finite element unknowns: " << size << endl;
+   }
+
+   // 7. Determine the list of true (i.e. parallel conforming) essential
+   //    boundary dofs. In this example, the boundary conditions are defined
+   //    by marking all the boundary attributes from the mesh as essential
+   //    (Dirichlet) and converting them to a list of true dofs.
+   Array<int> ess_tdof_list;
+   if (pmesh->bdr_attributes.Size())
+   {
+      Array<int> ess_bdr(pmesh->bdr_attributes.Max());
+      ess_bdr = 1;
+      fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+
+   // 8. Set up the parallel linear form b(.) which corresponds to the
+   //    right-hand side of the FEM linear system, which in this case is
+   //    (f,phi_i) where f is given by the function f_exact and phi_i are the
+   //    basis functions in the finite element fespace.
+   VectorFunctionCoefficient f(sdim, f_exact);
+   ParLinearForm *b = new ParLinearForm(fespace);
+   b->AddDomainIntegrator(new VectorFEDomainLFIntegrator(f));
+   b->Assemble();
+
+   // 9. Define the solution vector x as a parallel finite element grid function
+   //    corresponding to fespace. Initialize x by projecting the exact
+   //    solution. Note that only values from the boundary edges will be used
+   //    when eliminating the non-homogeneous boundary condition to modify the
+   //    r.h.s. vector b.
+   ParGridFunction x(fespace);
+   VectorFunctionCoefficient E(sdim, E_exact);
+   x.ProjectCoefficient(E);
+
+   // 10. Set up the parallel bilinear form corresponding to the EM diffusion
+   //     operator curl muinv curl + sigma I, by adding the curl-curl and the
+   //     mass domain integrators.
+   Coefficient *muinv = new ConstantCoefficient(1.0);
+   Coefficient *sigma = new ConstantCoefficient(-1.0);
+   ParBilinearForm *a = new ParBilinearForm(fespace);
+   a->AddDomainIntegrator(new CurlCurlIntegrator(*muinv));
+   a->AddDomainIntegrator(new VectorFEMassIntegrator(*sigma));
+
+   // 11. Assemble the parallel bilinear form and the corresponding linear
+   //     system, applying any necessary transformations such as: parallel
+   //     assembly, eliminating boundary conditions, applying conforming
+   //     constraints for non-conforming AMR, static condensation, etc.
+   if (static_cond) { a->EnableStaticCondensation(); }
+   a->Assemble();
+
+   HypreParMatrix A;
+   Vector B, X;
+   a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);
+
+   if (myid == 0)
+   {
+      cout << "Size of linear system: " << A.GetGlobalNumRows() << endl;
+   }
+
+   StopWatch chrono;
+   chrono.Clear();
+   chrono.Start();
+
+#ifdef MFEM_USE_STRUMPACK
+   if (use_strumpack)
+   {
+      Operator * Arow = new STRUMPACKRowLocMatrix(A);
+
+      STRUMPACKSolver * strumpack = new STRUMPACKSolver(argc, argv, MPI_COMM_WORLD);
+      strumpack->SetPrintFactorStatistics(true);
+      strumpack->SetPrintSolveStatistics(false);
+      strumpack->SetKrylovSolver(strumpack::KrylovSolver::DIRECT);
+      strumpack->SetReorderingStrategy(strumpack::ReorderingStrategy::METIS);
+      // strumpack->SetMC64Job(strumpack::MC64Job::NONE);
+      // strumpack->SetSymmetricPattern(true);
+      strumpack->SetOperator(*Arow);
+      strumpack->SetFromCommandLine();
+      //Solver * precond = strumpack;
+
+      strumpack->Mult(B, X);
+
+      delete strumpack;
+      delete Arow;
+   }
+   else
+#endif
+   {
+      // 12. Define and apply a parallel PCG solver for AX=B with the AMS
+      //     preconditioner from hypre.
+      ParFiniteElementSpace *prec_fespace =
+         (a->StaticCondensationIsEnabled() ? a->SCParFESpace() : fespace);
+      HypreSolver *ams = new HypreAMS(A, prec_fespace);
+      HyprePCG *pcg = new HyprePCG(A);
+      pcg->SetTol(1e-12);
+      pcg->SetMaxIter(500);
+      pcg->SetPrintLevel(2);
+      pcg->SetPreconditioner(*ams);
+      pcg->Mult(B, X);
+
+      delete pcg;
+      delete ams;
+   }
+
+   chrono.Stop();
+   cout << "Solver time " << chrono.RealTime() << endl;
+
+   // 13. Recover the parallel grid function corresponding to X. This is the
+   //     local finite element solution on each processor.
+   a->RecoverFEMSolution(X, *b, x);
+
+   // 14. Compute and print the L^2 norm of the error.
+   {
+      double err = x.ComputeL2Error(E);
+      if (myid == 0)
+      {
+         cout << "\n|| E_h - E ||_{L^2} = " << err << '\n' << endl;
+      }
+   }
+
+   // 15. Save the refined mesh and the solution in parallel. This output can
+   //     be viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
+   {
+      ostringstream mesh_name, sol_name;
+      mesh_name << "mesh." << setfill('0') << setw(6) << myid;
+      sol_name << "sol." << setfill('0') << setw(6) << myid;
+
+      ofstream mesh_ofs(mesh_name.str().c_str());
+      mesh_ofs.precision(8);
+      pmesh->Print(mesh_ofs);
+
+      ofstream sol_ofs(sol_name.str().c_str());
+      sol_ofs.precision(8);
+      x.Save(sol_ofs);
+   }
+
+   // 16. Send the solution by socket to a GLVis server.
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream sol_sock(vishost, visport);
+      sol_sock << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock.precision(8);
+      sol_sock << "solution\n" << *pmesh << x << flush;
+   }
+
+   // 17. Free the used memory.
+   delete a;
+   delete sigma;
+   delete muinv;
+   delete b;
+   delete fespace;
+   delete fec;
+   delete pmesh;
+
+   MPI_Finalize();
+
+   return 0;
+}
+
+
+void E_exact(const Vector &x, Vector &E)
+{
+   if (dim == 3)
+   {
+      E(0) = sin(kappa * x(1));
+      E(1) = sin(kappa * x(2));
+      E(2) = sin(kappa * x(0));
+   }
+   else
+   {
+      E(0) = sin(kappa * x(1));
+      E(1) = sin(kappa * x(0));
+      if (x.Size() == 3) { E(2) = 0.0; }
+   }
+}
+
+void f_exact(const Vector &x, Vector &f)
+{
+   if (dim == 3)
+   {
+      f(0) = (1. + kappa * kappa) * sin(kappa * x(1));
+      f(1) = (1. + kappa * kappa) * sin(kappa * x(2));
+      f(2) = (1. + kappa * kappa) * sin(kappa * x(0));
+   }
+   else
+   {
+      f(0) = (1. + kappa * kappa) * sin(kappa * x(1));
+      f(1) = (1. + kappa * kappa) * sin(kappa * x(0));
+      if (x.Size() == 3) { f(2) = 0.0; }
+   }
+}

examples/makefile

@@ -22,9 +22,9 @@ MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 
 SEQ_EXAMPLES = ex1 ex2 ex3 ex4 ex5 ex6 ex7 ex8 ex9 ex10 ex14 ex15 ex16 ex17\
-  ex18 ex19 ex20 ex22
+  ex18 ex19 ex20 ex21 ex22
 PAR_EXAMPLES = ex1p ex2p ex3p ex4p ex5p ex6p ex7p ex8p ex9p ex10p ex11p ex12p\
- ex13p ex14p ex15p ex16p ex17p ex18p ex19p ex20p ex22p
+ ex13p ex14p ex15p ex16p ex17p ex18p ex19p ex20p ex21p ex22p
 
 ifeq ($(MFEM_USE_MPI),NO)
    EXAMPLES = $(SEQ_EXAMPLES)
@@ -118,7 +118,7 @@ clean-build:
 clean-exec:
 	@rm -f refined.mesh displaced.mesh mesh.* ex5.mesh
 	@rm -rf Example5* Example9* Example15* Example16*
-	@rm -f sphere_refined.* sol.* sol_u.* sol_p.*
+	@rm -f sphere_refined.* sol.* sol_u.* sol_p.* sol_r.* sol_i.*
 	@rm -f ex9.mesh ex9-mesh.* ex9-init.* ex9-final.*
 	@rm -f deformed.* velocity.* elastic_energy.* mode_*
 	@rm -f ex16.mesh ex16-mesh.* ex16-init.* ex16-final.*

fem/complex_fem.cpp

@@ -0,0 +1,784 @@
+// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
+// reserved. See file COPYRIGHT for details.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability see http://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License (as published by the Free
+// Software Foundation) version 2.1 dated February 1999.
+
+#include "complex_fem.hpp"
+
+using namespace std;
+
+namespace mfem
+{
+
+ComplexGridFunction::ComplexGridFunction(FiniteElementSpace *fes)
+   : Vector(2*(fes->GetVSize()))
+{
+   gfr_ = new GridFunction(fes, &data[0]);
+   gfi_ = new GridFunction(fes, &data[fes->GetVSize()]);
+}
+
+void
+ComplexGridFunction::Update()
+{
+   FiniteElementSpace * fes = gfr_->FESpace();
+
+   int vsize = fes->GetVSize();
+
+   const Operator *T = fes->GetUpdateOperator();
+   if (T)
+   {
+      // Update the individual GridFunction objects.  This will allocate
+      // new data arrays for each GridFunction.
+      gfr_->Update();
+      gfi_->Update();
+
+      // Our data array now contains old data as well as being the wrong size
+      // so reallocate it.
+      this->SetSize(2 * vsize);
+
+      // Create temporary vectors which point to the new data array
+      Vector gf_r(&data[0], vsize);
+      Vector gf_i(&data[vsize], vsize);
+
+      // Copy the updated GridFunctions into the new data array
+      gf_r = *gfr_;
+      gf_i = *gfi_;
+
+      // Replace the individual data arrays with pointers into the new data array
+      gfr_->NewDataAndSize(&data[0], vsize);
+      gfi_->NewDataAndSize(&data[vsize], vsize);
+   }
+   else
+   {
+      // The existing data will not be transferred to the new GridFunctions
+      // so delete it a allocate a new array
+      this->SetSize(2 * vsize);
+
+      // Point the individual GridFunctions to the new data array
+      gfr_->NewDataAndSize(&data[0], vsize);
+      gfi_->NewDataAndSize(&data[vsize], vsize);
+
+      // These updates will only set the proper 'sequence' value within
+      // the individual GridFunction objects because their sizes are
+      // already correct
+      gfr_->Update();
+      gfi_->Update();
+   }
+}
+
+void
+ComplexGridFunction::ProjectCoefficient(Coefficient &real_coeff,
+                                        Coefficient &imag_coeff)
+{
+   gfr_->ProjectCoefficient(real_coeff);
+   gfi_->ProjectCoefficient(imag_coeff);
+}
+
+void
+ComplexGridFunction::ProjectCoefficient(VectorCoefficient &real_vcoeff,
+                                        VectorCoefficient &imag_vcoeff)
+{
+   gfr_->ProjectCoefficient(real_vcoeff);
+   gfi_->ProjectCoefficient(imag_vcoeff);
+}
+
+void
+ComplexGridFunction::ProjectBdrCoefficient(Coefficient &real_coeff,
+                                           Coefficient &imag_coeff,
+                                           Array<int> &attr)
+{
+   gfr_->ProjectBdrCoefficient(real_coeff, attr);
+   gfi_->ProjectBdrCoefficient(imag_coeff, attr);
+}
+
+void
+ComplexGridFunction::ProjectBdrCoefficientNormal(VectorCoefficient &real_vcoeff,
+                                                 VectorCoefficient &imag_vcoeff,
+                                                 Array<int> &attr)
+{
+   gfr_->ProjectBdrCoefficientNormal(real_vcoeff, attr);
+   gfi_->ProjectBdrCoefficientNormal(imag_vcoeff, attr);
+}
+
+void
+ComplexGridFunction::ProjectBdrCoefficientTangent(VectorCoefficient
+                                                  &real_vcoeff,
+                                                  VectorCoefficient
+                                                  &imag_vcoeff,
+                                                  Array<int> &attr)
+{
+   gfr_->ProjectBdrCoefficientTangent(real_vcoeff, attr);
+   gfi_->ProjectBdrCoefficientTangent(imag_vcoeff, attr);
+}
+
+
+ComplexLinearForm::ComplexLinearForm(FiniteElementSpace *f,
+                                     ComplexOperator::Convention convention)
+   : Vector(2*(f->GetVSize())),
+     conv_(convention)
+{
+   lfr_ = new LinearForm(f, &data[0]);
+   lfi_ = new LinearForm(f, &data[f->GetVSize()]);
+}
+
+ComplexLinearForm::~ComplexLinearForm()
+{
+   delete lfr_;
+   delete lfi_;
+}
+
+void
+ComplexLinearForm::AddDomainIntegrator(LinearFormIntegrator *lfi_real,
+                                       LinearFormIntegrator *lfi_imag)
+{
+   if ( lfi_real ) { lfr_->AddDomainIntegrator(lfi_real); }
+   if ( lfi_imag ) { lfi_->AddDomainIntegrator(lfi_imag); }
+}
+
+void
+ComplexLinearForm::Update()
+{
+   FiniteElementSpace *fes = lfr_->FESpace();
+
+   this->Update(fes);
+}
+
+void
+ComplexLinearForm::Update(FiniteElementSpace *fes)
+{
+   int vsize = fes->GetVSize();
+   SetSize(2 * vsize);
+
+   Vector lfr(&data[0], vsize);
+   Vector lfi(&data[vsize], vsize);
+
+   lfr_->Update(fes, lfr, 0);
+   lfi_->Update(fes, lfi, 0);
+}
+
+void
+ComplexLinearForm::Assemble()
+{
+   lfr_->Assemble();
+   lfi_->Assemble();
+   if (conv_ == ComplexOperator::BLOCK_SYMMETRIC)
+   {
+      *lfi_ *= -1.0;
+   }
+}
+
+complex<double>
+ComplexLinearForm::operator()(const ComplexGridFunction &gf) const
+{
+   double s = (conv_ == ComplexOperator::HERMITIAN)?1.0:-1.0;
+   return complex<double>((*lfr_)(gf.real()) - s * (*lfi_)(gf.imag()),
+                          (*lfr_)(gf.imag()) + s * (*lfi_)(gf.real()));
+}
+
+
+SesquilinearForm::SesquilinearForm(FiniteElementSpace *f,
+                                   ComplexOperator::Convention convention)
+   : conv_(convention),
+     blfr_(new BilinearForm(f)),
+     blfi_(new BilinearForm(f))
+{}
+
+SesquilinearForm::~SesquilinearForm()
+{
+   delete blfr_;
+   delete blfi_;
+}
+
+void SesquilinearForm::AddDomainIntegrator(BilinearFormIntegrator *bfi_real,
+                                           BilinearFormIntegrator *bfi_imag)
+{
+   if (bfi_real) { blfr_->AddDomainIntegrator(bfi_real); }
+   if (bfi_imag) { blfi_->AddDomainIntegrator(bfi_imag); }
+}
+
+void
+SesquilinearForm::AddBoundaryIntegrator(BilinearFormIntegrator *bfi_real,
+                                        BilinearFormIntegrator *bfi_imag)
+{
+   if (bfi_real) { blfr_->AddBoundaryIntegrator(bfi_real); }
+   if (bfi_imag) { blfi_->AddBoundaryIntegrator(bfi_imag); }
+}
+
+void
+SesquilinearForm::AddBoundaryIntegrator(BilinearFormIntegrator *bfi_real,
+                                        BilinearFormIntegrator *bfi_imag,
+                                        Array<int> & bdr_marker)
+{
+   if (bfi_real) { blfr_->AddBoundaryIntegrator(bfi_real, bdr_marker); }
+   if (bfi_imag) { blfi_->AddBoundaryIntegrator(bfi_imag, bdr_marker); }
+}
+
+void
+SesquilinearForm::Assemble(int skip_zeros)
+{
+   blfr_->Assemble(skip_zeros);
+   blfi_->Assemble(skip_zeros);
+}
+
+void
+SesquilinearForm::Finalize(int skip_zeros)
+{
+   blfr_->Finalize(skip_zeros);
+   blfi_->Finalize(skip_zeros);
+}
+
+ComplexSparseMatrix *
+SesquilinearForm::AssembleCompSpMat()
+{
+   return new ComplexSparseMatrix(&blfr_->SpMat(),
+                                  &blfi_->SpMat(),
+                                  false, false, conv_);
+
+}
+
+void
+SesquilinearForm::FormLinearSystem(const Array<int> &ess_tdof_list,
+                                   Vector &x, Vector &b,
+                                   OperatorHandle &A,
+                                   Vector &X, Vector &B,
+                                   int ci)
+{
+   FiniteElementSpace * fes = blfr_->FESpace();
+
+   int vsize  = fes->GetVSize();
+   // int tvsize = pfes->GetTrueVSize();
+
+   double s = (conv_ == ComplexOperator::HERMITIAN)?1.0:-1.0;
+
+   // Allocate temporary vectors
+   Vector b_0(vsize);  b_0 = 0.0;
+   // Vector B_0(tvsize); B_0 = 0.0;
+
+   // Extract the real and imaginary parts of the input vectors
+   MFEM_ASSERT(x.Size() == 2 * vsize, "Input GridFunction of incorrect size!");
+   Vector x_r(x.GetData(), vsize);
+   Vector x_i(&(x.GetData())[vsize], vsize);
+
+   MFEM_ASSERT(b.Size() == 2 * vsize, "Input LinearForm of incorrect size!");
+   Vector b_r(b.GetData(), vsize);
+   Vector b_i(&(b.GetData())[vsize], vsize);
+   b_i *= s;
+   /*
+   X.SetSize(2 * tvsize);
+   Vector X_r(X.GetData(), tvsize);
+   Vector X_i(&(X.GetData())[tvsize], tvsize);
+
+   B.SetSize(2 * tvsize);
+   Vector B_r(B.GetData(), tvsize);
+   Vector B_i(&(B.GetData())[tvsize], tvsize);
+   */
+   SparseMatrix * A_r = new SparseMatrix;
+   SparseMatrix * A_i = new SparseMatrix;
+   Vector X_0, B_0;
+
+   b_0 = b_r;
+   blfr_->FormLinearSystem(ess_tdof_list, x_r, b_r, *A_r, X_0, B_0, ci);
+
+   int tvsize = B_0.Size();
+   X.SetSize(2 * tvsize);
+   B.SetSize(2 * tvsize);
+   Vector X_r(X.GetData(), tvsize);
+   Vector X_i(&(X.GetData())[tvsize], tvsize);
+   Vector B_r(B.GetData(), tvsize);
+   Vector B_i(&(B.GetData())[tvsize], tvsize);
+   X_r = X_0; B_r = B_0;
+
+   b_0 = 0.0;
+   blfi_->FormLinearSystem(ess_tdof_list, x_i, b_0, *A_i, X_0, B_0, false);
+   B_r -= B_0;
+
+   b_0 = b_i;
+   blfr_->FormLinearSystem(ess_tdof_list, x_i, b_0, *A_r, X_0, B_0, ci);
+   X_i = X_0; B_i = B_0;
+
+   b_0 = 0.0;
+   blfi_->FormLinearSystem(ess_tdof_list, x_r, b_0, *A_i, X_0, B_0, false);
+   B_i += B_0;
+
+   B_i *= s;
+   b_i *= s;
+
+   // A = A_r + i A_i
+   A.Clear();
+   ComplexSparseMatrix * A_sp =
+      new ComplexSparseMatrix(A_r, A_i, true, true, conv_);
+   A.Reset<ComplexSparseMatrix>(A_sp, true);
+}
+
+void
+SesquilinearForm::RecoverFEMSolution(const Vector &X, const Vector &b,
+                                     Vector &x)
+{
+   FiniteElementSpace * fes = blfr_->FESpace();
+
+   const SparseMatrix *P = fes->GetConformingProlongation();
+
+   int vsize  = fes->GetVSize();
+   int tvsize = X.Size() / 2;
+
+   Vector X_r(X.GetData(), tvsize);
+   Vector X_i(&(X.GetData())[tvsize], tvsize);
+
+   Vector x_r(x.GetData(), vsize);
+   Vector x_i(&(x.GetData())[vsize], vsize);
+
+   if (!P)
+   {
+      x = X;
+   }
+   else
+   {
+      // Apply conforming prolongation
+      P->Mult(X_r, x_r);
+      P->Mult(X_i, x_i);
+   }
+}
+
+void
+SesquilinearForm::Update(FiniteElementSpace *nfes)
+{
+   if ( blfr_ ) { blfr_->Update(nfes); }
+   if ( blfi_ ) { blfi_->Update(nfes); }
+}
+
+
+#ifdef MFEM_USE_MPI
+
+ParComplexGridFunction::ParComplexGridFunction(ParFiniteElementSpace *pfes)
+   : Vector(2*(pfes->GetVSize()))
+{
+   pgfr_ = new ParGridFunction(pfes, &data[0]);
+   pgfi_ = new ParGridFunction(pfes, &data[pfes->GetVSize()]);
+}
+
+void
+ParComplexGridFunction::Update()
+{
+   ParFiniteElementSpace * pfes = pgfr_->ParFESpace();
+
+   int vsize = pfes->GetVSize();
+
+   const Operator *T = pfes->GetUpdateOperator();
+   if (T)
+   {
+      // Update the individual GridFunction objects.  This will allocate
+      // new data arrays for each GridFunction.
+      pgfr_->Update();
+      pgfi_->Update();
+
+      // Our data array now contains old data as well as being the wrong size
+      // so reallocate it.
+      this->SetSize(2 * vsize);
+
+      // Create temporary vectors which point to the new data array
+      Vector gf_r(&data[0], vsize);
+      Vector gf_i(&data[vsize], vsize);
+
+      // Copy the updated GridFunctions into the new data array
+      gf_r = *pgfr_;
+      gf_i = *pgfi_;
+
+      // Replace the individual data arrays with pointers into the new data array
+      pgfr_->NewDataAndSize(&data[0], vsize);
+      pgfi_->NewDataAndSize(&data[vsize], vsize);
+   }
+   else
+   {
+      // The existing data will not be transferred to the new GridFunctions
+      // so delete it a allocate a new array
+      this->SetSize(2 * vsize);
+
+      // Point the individual GridFunctions to the new data array
+      pgfr_->NewDataAndSize(&data[0], vsize);
+      pgfi_->NewDataAndSize(&data[vsize], vsize);
+
+      // These updates will only set the proper 'sequence' value within
+      // the individual GridFunction objects because their sizes are
+      // already correct
+      pgfr_->Update();
+      pgfi_->Update();
+   }
+}
+
+void
+ParComplexGridFunction::ProjectCoefficient(Coefficient &real_coeff,
+                                           Coefficient &imag_coeff)
+{
+   pgfr_->ProjectCoefficient(real_coeff);
+   pgfi_->ProjectCoefficient(imag_coeff);
+}
+
+void
+ParComplexGridFunction::ProjectCoefficient(VectorCoefficient &real_vcoeff,
+                                           VectorCoefficient &imag_vcoeff)
+{
+   pgfr_->ProjectCoefficient(real_vcoeff);
+   pgfi_->ProjectCoefficient(imag_vcoeff);
+}
+
+void
+ParComplexGridFunction::ProjectBdrCoefficient(Coefficient &real_coeff,
+                                              Coefficient &imag_coeff,
+                                              Array<int> &attr)
+{
+   pgfr_->ProjectBdrCoefficient(real_coeff, attr);
+   pgfi_->ProjectBdrCoefficient(imag_coeff, attr);
+}
+
+void
+ParComplexGridFunction::ProjectBdrCoefficientNormal(VectorCoefficient
+                                                    &real_vcoeff,
+                                                    VectorCoefficient
+                                                    &imag_vcoeff,
+                                                    Array<int> &attr)
+{
+   pgfr_->ProjectBdrCoefficientNormal(real_vcoeff, attr);
+   pgfi_->ProjectBdrCoefficientNormal(imag_vcoeff, attr);
+}
+
+void
+ParComplexGridFunction::ProjectBdrCoefficientTangent(VectorCoefficient
+                                                     &real_vcoeff,
+                                                     VectorCoefficient
+                                                     &imag_vcoeff,
+                                                     Array<int> &attr)
+{
+   pgfr_->ProjectBdrCoefficientTangent(real_vcoeff, attr);
+   pgfi_->ProjectBdrCoefficientTangent(imag_vcoeff, attr);
+}
+
+void
+ParComplexGridFunction::Distribute(const Vector *tv)
+{
+   ParFiniteElementSpace * pfes = pgfr_->ParFESpace();
+   HYPRE_Int size = pfes->GetTrueVSize();
+
+   double * tvd = tv->GetData();
+   Vector tvr(tvd, size);
+   Vector tvi(&tvd[size], size);
+
+   pgfr_->Distribute(tvr);
+   pgfi_->Distribute(tvi);
+}
+
+void
+ParComplexGridFunction::ParallelProject(Vector &tv) const
+{
+   ParFiniteElementSpace * pfes = pgfr_->ParFESpace();
+   HYPRE_Int size = pfes->GetTrueVSize();
+
+   double * tvd = tv.GetData();
+   Vector tvr(tvd, size);
+   Vector tvi(&tvd[size], size);
+
+   pgfr_->ParallelProject(tvr);
+   pgfi_->ParallelProject(tvi);
+}
+
+
+ParComplexLinearForm::ParComplexLinearForm(ParFiniteElementSpace *pfes,
+                                           ComplexOperator::Convention
+                                           convention)
+   : Vector(2*(pfes->GetVSize())),
+     conv_(convention)
+{
+   plfr_ = new ParLinearForm(pfes, &data[0]);
+   plfi_ = new ParLinearForm(pfes, &data[pfes->GetVSize()]);
+
+   HYPRE_Int * tdof_offsets = pfes->GetTrueDofOffsets();
+
+   int n = (HYPRE_AssumedPartitionCheck()) ? 2 : pfes->GetNRanks();
+   tdof_offsets_ = new HYPRE_Int[n+1];
+
+   for (int i=0; i<=n; i++)
+   {
+      tdof_offsets_[i] = 2 * tdof_offsets[i];
+   }
+}
+
+ParComplexLinearForm::~ParComplexLinearForm()
+{
+   delete plfr_;
+   delete plfi_;
+   delete [] tdof_offsets_;
+}
+
+void
+ParComplexLinearForm::AddDomainIntegrator(LinearFormIntegrator *lfi_real,
+                                          LinearFormIntegrator *lfi_imag)
+{
+   if ( lfi_real ) { plfr_->AddDomainIntegrator(lfi_real); }
+   if ( lfi_imag ) { plfi_->AddDomainIntegrator(lfi_imag); }
+}
+
+void
+ParComplexLinearForm::Update(ParFiniteElementSpace *pf)
+{
+   ParFiniteElementSpace *pfes = (pf!=NULL)?pf:plfr_->ParFESpace();
+   int vsize = pfes->GetVSize();
+   SetSize(2 * vsize);
+
+   Vector plfr(&data[0], vsize);
+   Vector plfi(&data[vsize], vsize);
+
+   plfr_->Update(pfes, plfr, 0);
+   plfi_->Update(pfes, plfi, 0);
+}
+
+void
+ParComplexLinearForm::Assemble()
+{
+   plfr_->Assemble();
+   plfi_->Assemble();
+   if (conv_ == ComplexOperator::BLOCK_SYMMETRIC)
+   {
+      *plfi_ *= -1.0;
+   }
+}
+
+void
+ParComplexLinearForm::ParallelAssemble(Vector &tv)
+{
+   HYPRE_Int size = plfr_->ParFESpace()->GetTrueVSize();
+
+   double * tvd = tv.GetData();
+   Vector tvr(tvd, size);
+   Vector tvi(&tvd[size], size);
+
+   plfr_->ParallelAssemble(tvr);
+   plfi_->ParallelAssemble(tvi);
+}
+
+HypreParVector *
+ParComplexLinearForm::ParallelAssemble()
+{
+   const ParFiniteElementSpace * pfes = plfr_->ParFESpace();
+
+   HypreParVector * tv = new HypreParVector(pfes->GetComm(),
+                                            2*(pfes->GlobalTrueVSize()),
+                                            tdof_offsets_);
+
+   HYPRE_Int size = pfes->GetTrueVSize();
+
+   double * tvd = tv->GetData();
+   Vector tvr(tvd, size);
+   Vector tvi(&tvd[size], size);
+
+   plfr_->ParallelAssemble(tvr);
+   plfi_->ParallelAssemble(tvi);
+
+   return tv;
+}
+
+complex<double>
+ParComplexLinearForm::operator()(const ParComplexGridFunction &gf) const
+{
+   double s = (conv_ == ComplexOperator::HERMITIAN)?1.0:-1.0;
+   return complex<double>((*plfr_)(gf.real()) - s * (*plfi_)(gf.imag()),
+                          (*plfr_)(gf.imag()) + s * (*plfi_)(gf.real()));
+}
+
+
+ParSesquilinearForm::ParSesquilinearForm(ParFiniteElementSpace *pf,
+                                         ComplexOperator::Convention
+                                         convention)
+   : conv_(convention),
+     pblfr_(new ParBilinearForm(pf)),
+     pblfi_(new ParBilinearForm(pf))
+{}
+
+ParSesquilinearForm::~ParSesquilinearForm()
+{
+   delete pblfr_;
+   delete pblfi_;
+}
+
+void ParSesquilinearForm::AddDomainIntegrator(BilinearFormIntegrator *bfi_real,
+                                              BilinearFormIntegrator *bfi_imag)
+{
+   if (bfi_real) { pblfr_->AddDomainIntegrator(bfi_real); }
+   if (bfi_imag) { pblfi_->AddDomainIntegrator(bfi_imag); }
+}
+
+void
+ParSesquilinearForm::AddBoundaryIntegrator(BilinearFormIntegrator *bfi_real,
+                                           BilinearFormIntegrator *bfi_imag)
+{
+   if (bfi_real) { pblfr_->AddBoundaryIntegrator(bfi_real); }
+   if (bfi_imag) { pblfi_->AddBoundaryIntegrator(bfi_imag); }
+}
+
+void
+ParSesquilinearForm::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
+ParSesquilinearForm::Assemble(int skip_zeros)
+{
+   pblfr_->Assemble(skip_zeros);
+   pblfi_->Assemble(skip_zeros);
+}
+
+void
+ParSesquilinearForm::Finalize(int skip_zeros)
+{
+   pblfr_->Finalize(skip_zeros);
+   pblfi_->Finalize(skip_zeros);
+}
+
+ComplexHypreParMatrix *
+ParSesquilinearForm::ParallelAssemble()
+{
+   return new ComplexHypreParMatrix(pblfr_->ParallelAssemble(),
+                                    pblfi_->ParallelAssemble(),
+                                    true, true, conv_);
+
+}
+
+void
+ParSesquilinearForm::FormLinearSystem(const Array<int> &ess_tdof_list,
+                                      Vector &x, Vector &b,
+                                      OperatorHandle &A,
+                                      Vector &X, Vector &B,
+                                      int ci)
+{
+   ParFiniteElementSpace * pfes = pblfr_->ParFESpace();
+
+   int tvs = pfes->TrueVSize();
+   cout << "TrueVSize returns " << tvs << endl;
+   cout << "GetVSize returns " << pfes->GetVSize() << endl;
+
+   int vsize = x.Size() / 2;
+   // int vsize  = pfes->GetVSize();
+   // int tvsize = pfes->GetTrueVSize();
+
+   cout << "x.Size/2 returns " << vsize << endl;
+
+   double s = (conv_ == ComplexOperator::HERMITIAN)?1.0:-1.0;
+
+   // Allocate temporary vectors
+   Vector b_0(vsize);  b_0 = 0.0;
+   // Vector B_0(tvsize); B_0 = 0.0;
+
+   // Extract the real and imaginary parts of the input vectors
+   // MFEM_ASSERT(x.Size() == 2 * vsize, "Input GridFunction of incorrect size!");
+   Vector x_r(x.GetData(), vsize);
+   Vector x_i(&(x.GetData())[vsize], vsize);
+
+   MFEM_ASSERT(b.Size() == 2 * vsize, "Input LinearForm of incorrect size!");
+   Vector b_r(b.GetData(), vsize);
+   Vector b_i(&(b.GetData())[vsize], vsize);
+   b_i *= s;
+   /*
+   X.SetSize(2 * tvsize);
+   Vector X_r(X.GetData(), tvsize);
+   Vector X_i(&(X.GetData())[tvsize], tvsize);
+
+   B.SetSize(2 * tvsize);
+   Vector B_r(B.GetData(), tvsize);
+   Vector B_i(&(B.GetData())[tvsize], tvsize);
+   */
+   OperatorHandle A_r, A_i;
+   Vector X_0, B_0;
+   cout << "pblfr fls 1" << endl << flush;
+   b_0 = b_r;
+   pblfr_->FormLinearSystem(ess_tdof_list, x_r, b_0, A_r, X_0, B_0, ci);
+
+   int tvsize = B_0.Size();
+   X.SetSize(2 * tvsize);
+   B.SetSize(2 * tvsize);
+   Vector X_r(X.GetData(), tvsize);
+   Vector X_i(&(X.GetData())[tvsize], tvsize);
+   Vector B_r(B.GetData(), tvsize);
+   Vector B_i(&(B.GetData())[tvsize], tvsize);
+   X_r = X_0; B_r = B_0;
+   cout << "pblfi fls 1" << endl << flush;
+   b_0 = 0.0;
+   pblfi_->FormLinearSystem(ess_tdof_list, x_i, b_0, A_i, X_0, B_0, false);
+   B_r -= B_0;
+   cout << "pblfr fls 2" << endl << flush;
+   b_0 = b_i;
+   pblfr_->FormLinearSystem(ess_tdof_list, x_i, b_0, A_r, X_0, B_0, ci);
+   X_i = X_0; B_i = B_0;
+   cout << "pblfi fls 2" << endl << flush;
+   b_0 = 0.0;
+   pblfi_->FormLinearSystem(ess_tdof_list, x_r, b_0, A_i, X_0, B_0, false);
+   B_i += B_0;
+
+   B_i *= s;
+   b_i *= s;
+
+   // 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);
+   }
+}
+
+void
+ParSesquilinearForm::RecoverFEMSolution(const Vector &X, const Vector &b,
+                                        Vector &x)
+{
+   ParFiniteElementSpace * pfes = pblfr_->ParFESpace();
+
+   const Operator &P = *pfes->GetProlongationMatrix();
+
+   int vsize  = pfes->GetVSize();
+   int tvsize = X.Size() / 2;
+
+   Vector X_r(X.GetData(), tvsize);
+   Vector X_i(&(X.GetData())[tvsize], tvsize);
+
+   Vector x_r(x.GetData(), vsize);
+   Vector x_i(&(x.GetData())[vsize], vsize);
+
+   // Apply conforming prolongation
+   P.Mult(X_r, x_r);
+   P.Mult(X_i, x_i);
+}
+
+void
+ParSesquilinearForm::Update(FiniteElementSpace *nfes)
+{
+   if ( pblfr_ ) { pblfr_->Update(nfes); }
+   if ( pblfi_ ) { pblfi_->Update(nfes); }
+}
+
+
+#endif // MFEM_USE_MPI
+
+}

fem/complex_fem.hpp

@@ -0,0 +1,356 @@
+// 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_COMPLEX_FEM
+#define MFEM_COMPLEX_FEM
+
+#include "../linalg/complex_operator.hpp"
+#include "gridfunc.hpp"
+#include "linearform.hpp"
+#include "bilinearform.hpp"
+#ifdef MFEM_USE_MPI
+#include "pgridfunc.hpp"
+#include "plinearform.hpp"
+#include "pbilinearform.hpp"
+#endif
+#include <complex>
+
+namespace mfem
+{
+
+/// Class for complex-valued grid function - Vector with associated FE space.
+class ComplexGridFunction : public Vector
+{
+private:
+
+   GridFunction * gfr_;
+   GridFunction * gfi_;
+
+protected:
+   void Destroy() { delete gfr_; delete gfi_; }
+
+public:
+
+   /* @brief Construct a ComplexGridFunction associated with the
+      FiniteElementSpace @a *f. */
+   ComplexGridFunction(FiniteElementSpace *f);
+
+   void Update();
+
+   /// Assign constant values to the ComplexGridFunction data.
+   ComplexGridFunction &operator=(const std::complex<double> & value)
+   { *gfr_ = value.real(); *gfi_ = value.imag(); return *this; }
+
+   virtual void ProjectCoefficient(Coefficient &real_coeff,
+                                   Coefficient &imag_coeff);
+   virtual void ProjectCoefficient(VectorCoefficient &real_vcoeff,
+                                   VectorCoefficient &imag_vcoeff);
+
+   virtual void ProjectBdrCoefficient(Coefficient &real_coeff,
+                                      Coefficient &imag_coeff,
+                                      Array<int> &attr);
+   virtual void ProjectBdrCoefficientNormal(VectorCoefficient &real_coeff,
+                                            VectorCoefficient &imag_coeff,
+                                            Array<int> &attr);
+   virtual void ProjectBdrCoefficientTangent(VectorCoefficient &real_coeff,
+                                             VectorCoefficient &imag_coeff,
+                                             Array<int> &attr);
+
+   FiniteElementSpace *FESpace() { return gfr_->FESpace(); }
+   const FiniteElementSpace *FESpace() const { return gfr_->FESpace(); }
+
+   GridFunction & real() { return *gfr_; }
+   GridFunction & imag() { return *gfi_; }
+   const GridFunction & real() const { return *gfr_; }
+   const GridFunction & imag() const { return *gfi_; }
+
+   /// Destroys grid function.
+   virtual ~ComplexGridFunction() { Destroy(); }
+
+};
+
+class ComplexLinearForm : public Vector
+{
+private:
+   ComplexOperator::Convention conv_;
+
+protected:
+   LinearForm * lfr_;
+   LinearForm * lfi_;
+
+   // HYPRE_Int * tdof_offsets_;
+
+public:
+
+   ComplexLinearForm(FiniteElementSpace *fes,
+                     ComplexOperator::Convention
+                     convention = ComplexOperator::HERMITIAN);
+
+   virtual ~ComplexLinearForm();
+
+   /// Adds new Domain Integrator.
+   void AddDomainIntegrator(LinearFormIntegrator *lfi_real,
+                            LinearFormIntegrator *lfi_imag);
+
+   FiniteElementSpace *FESpace() const { return lfr_->FESpace(); }
+
+   LinearForm & real() { return *lfr_; }
+   LinearForm & imag() { return *lfi_; }
+   const LinearForm & real() const { return *lfr_; }
+   const LinearForm & imag() const { return *lfi_; }
+
+   void Update();
+   void Update(FiniteElementSpace *f);
+
+   /// Assembles the linear form i.e. sums over all domain/bdr integrators.
+   void Assemble();
+
+   std::complex<double> operator()(const ComplexGridFunction &gf) const;
+
+};
+
+// Class for sesquilinear form
+class SesquilinearForm
+{
+private:
+   ComplexOperator::Convention conv_;
+
+   //protected:
+   BilinearForm *blfr_;
+   BilinearForm *blfi_;
+
+public:
+   SesquilinearForm(FiniteElementSpace *fes,
+                    ComplexOperator::Convention
+                    convention = ComplexOperator::HERMITIAN);
+
+   ComplexOperator::Convention GetConvention() const { return conv_; }
+   void SetConvention(const ComplexOperator::Convention &
+                      convention) { conv_  = convention; }
+
+   BilinearForm & real() { return *blfr_; }
+   BilinearForm & imag() { return *blfi_; }
+   const BilinearForm & real() const { return *blfr_; }
+   const BilinearForm & imag() const { return *blfi_; }
+
+   /// Adds new Domain Integrator.
+   void AddDomainIntegrator(BilinearFormIntegrator *bfi_real,
+                            BilinearFormIntegrator *bfi_imag);
+
+   /// Adds new Boundary Integrator.
+   void AddBoundaryIntegrator(BilinearFormIntegrator *bfi_real,
+                              BilinearFormIntegrator *bfi_imag);
+
+   /// Adds new Boundary Integrator, restricted to specific boundary attributes.
+   void AddBoundaryIntegrator(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^t A P.
+   /** The returned matrix has to be deleted by the caller. */
+   ComplexSparseMatrix *AssembleCompSpMat();
+
+   /// Return the parallel FE space associated with the ParBilinearForm.
+   FiniteElementSpace *FESpace() const { return blfr_->FESpace(); }
+
+   void FormLinearSystem(const Array<int> &ess_tdof_list, Vector &x, Vector &b,
+                         OperatorHandle &A, Vector &X, Vector &B,
+                         int copy_interior = 0);
+
+   /** Call this method after solving a linear system constructed using the
+       FormLinearSystem method to recover the solution as a ParGridFunction-size
+       vector in x. Use the same arguments as in the FormLinearSystem call. */
+   virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x);
+
+   virtual void Update(FiniteElementSpace *nfes = NULL);
+
+   virtual ~SesquilinearForm();
+};
+
+#ifdef MFEM_USE_MPI
+
+/// Class for complex-valued grid function - Vector with associated FE space.
+class ParComplexGridFunction : public Vector
+{
+private:
+
+   ParGridFunction * pgfr_;
+   ParGridFunction * pgfi_;
+
+protected:
+   void Destroy() { delete pgfr_; delete pgfi_; }
+
+public:
+
+   /* @brief Construct a ParComplexGridFunction associated with the
+      ParFiniteElementSpace @a *f. */
+   ParComplexGridFunction(ParFiniteElementSpace *pf);
+
+   void Update();
+
+   /// Assign constant values to the ParComplexGridFunction data.
+   ParComplexGridFunction &operator=(const std::complex<double> & value)
+   { *pgfr_ = value.real(); *pgfi_ = value.imag(); return *this; }
+
+   virtual void ProjectCoefficient(Coefficient &real_coeff,
+                                   Coefficient &imag_coeff);
+   virtual void ProjectCoefficient(VectorCoefficient &real_vcoeff,
+                                   VectorCoefficient &imag_vcoeff);
+
+   virtual void ProjectBdrCoefficient(Coefficient &real_coeff,
+                                      Coefficient &imag_coeff,
+                                      Array<int> &attr);
+   virtual void ProjectBdrCoefficientNormal(VectorCoefficient &real_coeff,
+                                            VectorCoefficient &imag_coeff,
+                                            Array<int> &attr);
+   virtual void ProjectBdrCoefficientTangent(VectorCoefficient &real_coeff,
+                                             VectorCoefficient &imag_coeff,
+                                             Array<int> &attr);
+
+   void Distribute(const Vector *tv);
+   void Distribute(const Vector &tv) { Distribute(&tv); }
+
+   /// Returns the vector restricted to the true dofs.
+   void ParallelProject(Vector &tv) const;
+
+   FiniteElementSpace *FESpace() { return pgfr_->FESpace(); }
+   const FiniteElementSpace *FESpace() const { return pgfr_->FESpace(); }
+
+   ParGridFunction & real() { return *pgfr_; }
+   ParGridFunction & imag() { return *pgfi_; }
+   const ParGridFunction & real() const { return *pgfr_; }
+   const ParGridFunction & imag() const { return *pgfi_; }
+
+   /// Destroys grid function.
+   virtual ~ParComplexGridFunction() { Destroy(); }
+
+};
+
+class ParComplexLinearForm : public Vector
+{
+private:
+   ComplexOperator::Convention conv_;
+
+protected:
+   ParLinearForm * plfr_;
+   ParLinearForm * plfi_;
+
+   HYPRE_Int * tdof_offsets_;
+
+public:
+
+   ParComplexLinearForm(ParFiniteElementSpace *pf,
+                        ComplexOperator::Convention
+                        convention = ComplexOperator::HERMITIAN);
+
+   virtual ~ParComplexLinearForm();
+
+   /// Adds new Domain Integrator.
+   void AddDomainIntegrator(LinearFormIntegrator *lfi_real,
+                            LinearFormIntegrator *lfi_imag);
+
+   ParFiniteElementSpace *ParFESpace() const { return plfr_->ParFESpace(); }
+
+   ParLinearForm & real() { return *plfr_; }
+   ParLinearForm & imag() { return *plfi_; }
+   const ParLinearForm & real() const { return *plfr_; }
+   const ParLinearForm & imag() const { return *plfi_; }
+
+   void Update(ParFiniteElementSpace *pf = NULL);
+
+   /// Assembles the linear form i.e. sums over all domain/bdr integrators.
+   void Assemble();
+
+   /// Assemble the vector on the true dofs, i.e. P^t v.
+   void ParallelAssemble(Vector &tv);
+
+   /// Returns the vector assembled on the true dofs, i.e. P^t v.
+   HypreParVector *ParallelAssemble();
+
+   std::complex<double> operator()(const ParComplexGridFunction &gf) const;
+
+};
+
+// Class for parallel sesquilinear form
+class ParSesquilinearForm
+{
+private:
+   ComplexOperator::Convention conv_;
+
+   //protected:
+   ParBilinearForm *pblfr_;
+   ParBilinearForm *pblfi_;
+
+public:
+   ParSesquilinearForm(ParFiniteElementSpace *pf,
+                       ComplexOperator::Convention
+                       convention = ComplexOperator::HERMITIAN);
+
+   ComplexOperator::Convention GetConvention() const { return conv_; }
+   void SetConvention(const ComplexOperator::Convention &
+                      convention) { conv_  = convention; }
+
+   ParBilinearForm & real() { return *pblfr_; }
+   ParBilinearForm & imag() { return *pblfi_; }
+   const ParBilinearForm & real() const { return *pblfr_; }
+   const ParBilinearForm & imag() const { return *pblfi_; }
+
+   /// Adds new Domain Integrator.
+   void AddDomainIntegrator(BilinearFormIntegrator *bfi_real,
+                            BilinearFormIntegrator *bfi_imag);
+
+   /// Adds new Boundary Integrator.
+   void AddBoundaryIntegrator(BilinearFormIntegrator *bfi_real,
+                              BilinearFormIntegrator *bfi_imag);
+
+   /// Adds new Boundary Integrator, restricted to specific boundary attributes.
+   void AddBoundaryIntegrator(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^t A P.
+   /** The returned matrix has to be deleted by the caller. */
+   ComplexHypreParMatrix *ParallelAssemble();
+
+   /// Return the parallel FE space associated with the ParBilinearForm.
+   ParFiniteElementSpace *ParFESpace() const { return pblfr_->ParFESpace(); }
+
+   void FormLinearSystem(const Array<int> &ess_tdof_list, Vector &x, Vector &b,
+                         OperatorHandle &A, Vector &X, Vector &B,
+                         int copy_interior = 0);
+
+   /** Call this method after solving a linear system constructed using the
+       FormLinearSystem method to recover the solution as a ParGridFunction-size
+       vector in x. Use the same arguments as in the FormLinearSystem call. */
+   virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x);
+
+   virtual void Update(FiniteElementSpace *nfes = NULL);
+
+   virtual ~ParSesquilinearForm();
+};
+
+#endif // MFEM_USE_MPI
+
+}
+
+#endif // MFEM_COMPLEX_FEM

fem/fem.hpp

@@ -18,6 +18,7 @@
 #include "fe_coll.hpp"
 #include "eltrans.hpp"
 #include "coefficient.hpp"
+#include "complex_fem.hpp"
 #include "lininteg.hpp"
 #include "nonlininteg.hpp"
 #include "bilininteg.hpp"

fem/linearform.hpp

@@ -81,6 +81,15 @@ public:
        Update(FiniteElementSpace *, Vector &, int). */
    LinearForm() { fes = NULL; extern_lfs = 0; }
 
+   /// Construct a LinearForm using previously allocated array @a data.
+   /** The LinearForm does not assume ownership of @a data which is assumed to
+       be of size at least `f->GetVSize()`. Similar to the Vector constructor
+       for externally allocated array, the pointer @a data can be NULL. The data
+       array can be replaced later using the method SetData().
+    */
+   LinearForm(FiniteElementSpace *f, double *data) : Vector(data, f->GetVSize())
+   { fes = f; }
+
    /// Copy assignment. Only the data of the base class Vector is copied.
    /** It is assumed that this object and @a rhs use FiniteElementSpace%s that
        have the same size.

fem/plinearform.hpp

@@ -45,6 +45,16 @@ public:
    /** The pointer @a pf is not owned by the newly constructed object. */
    ParLinearForm(ParFiniteElementSpace *pf) : LinearForm(pf) { pfes = pf; }
 
+   /// Construct a ParLinearForm using previously allocated array @a data.
+   /** The ParLinearForm does not assume ownership of @a data which is assumed
+       to be of size at least `pf->GetVSize()`. Similar to the LinearForm and
+       Vector constructors for externally allocated array, the pointer @a data
+       can be NULL. The data array can be replaced later using the method
+       SetData().
+    */
+   ParLinearForm(ParFiniteElementSpace *pf, double *data) :
+      LinearForm(pf, data), pfes(pf) { }
+
    /** @brief Create a ParLinearForm on the ParFiniteElementSpace @a *pf, using
        the same integrators as the ParLinearForm @a *plf.
 

linalg/complex_operator.cpp

@@ -10,6 +10,8 @@
 // Software Foundation) version 2.1 dated February 1999.
 
 #include "complex_operator.hpp"
+#include <set>
+#include <map>
 
 namespace mfem
 {
@@ -39,6 +41,30 @@ ComplexOperator::~ComplexOperator()
    delete v_;
 }
 
+Operator & ComplexOperator::real()
+{
+   MFEM_ASSERT(Op_Real_, "ComplexOperator has no real part!");
+   return *Op_Real_;
+}
+
+Operator & ComplexOperator::imag()
+{
+   MFEM_ASSERT(Op_Imag_, "ComplexOperator has no imaginary part!");
+   return *Op_Imag_;
+}
+
+const Operator & ComplexOperator::real() const
+{
+   MFEM_ASSERT(Op_Real_, "ComplexOperator has no real part!");
+   return *Op_Real_;
+}
+
+const Operator & ComplexOperator::imag() const
+{
+   MFEM_ASSERT(Op_Imag_, "ComplexOperator has no imaginary part!");
+   return *Op_Imag_;
+}
+
 void ComplexOperator::Mult(const Vector &x, Vector &y) const
 {
    double * x_data = x.GetData();
@@ -120,6 +146,30 @@ void ComplexOperator::MultTranspose(const Vector &x_r, const Vector &x_i,
 }
 
 
+SparseMatrix & ComplexSparseMatrix::real()
+{
+   MFEM_ASSERT(Op_Real_, "ComplexSparseMatrix has no real part!");
+   return dynamic_cast<SparseMatrix &>(*Op_Real_);
+}
+
+SparseMatrix & ComplexSparseMatrix::imag()
+{
+   MFEM_ASSERT(Op_Imag_, "ComplexSparseMatrix has no imaginary part!");
+   return dynamic_cast<SparseMatrix &>(*Op_Imag_);
+}
+
+const SparseMatrix & ComplexSparseMatrix::real() const
+{
+   MFEM_ASSERT(Op_Real_, "ComplexSparseMatrix has no real part!");
+   return dynamic_cast<const SparseMatrix &>(*Op_Real_);
+}
+
+const SparseMatrix & ComplexSparseMatrix::imag() const
+{
+   MFEM_ASSERT(Op_Imag_, "ComplexSparseMatrix has no imaginary part!");
+   return dynamic_cast<const SparseMatrix &>(*Op_Imag_);
+}
+
 SparseMatrix * ComplexSparseMatrix::GetSystemMatrix() const
 {
    SparseMatrix * A_r = dynamic_cast<SparseMatrix*>(Op_Real_);
@@ -184,4 +234,358 @@ SparseMatrix * ComplexSparseMatrix::GetSystemMatrix() const
    return new SparseMatrix(I, J, D, this->Height(), this->Width());
 }
 
+#ifdef MFEM_USE_MPI
+
+ComplexHypreParMatrix::ComplexHypreParMatrix(HypreParMatrix * A_Real,
+                                             HypreParMatrix * A_Imag,
+                                             bool ownReal, bool ownImag,
+                                             Convention convention)
+   : ComplexOperator(A_Real, A_Imag, ownReal, ownImag, convention)
+{
+   comm_ = (A_Real) ? A_Real->GetComm() :
+           ((A_Imag) ? A_Imag->GetComm() : MPI_COMM_WORLD);
+
+   MPI_Comm_rank(comm_, &myid_);
+   MPI_Comm_size(comm_, &nranks_);
+}
+
+HypreParMatrix & ComplexHypreParMatrix::real()
+{
+   MFEM_ASSERT(Op_Real_, "ComplexHypreParMatrix has no real part!");
+   return dynamic_cast<HypreParMatrix &>(*Op_Real_);
+}
+
+HypreParMatrix & ComplexHypreParMatrix::imag()
+{
+   MFEM_ASSERT(Op_Imag_, "ComplexHypreParMatrix has no imaginary part!");
+   return dynamic_cast<HypreParMatrix &>(*Op_Imag_);
+}
+
+const HypreParMatrix & ComplexHypreParMatrix::real() const
+{
+   MFEM_ASSERT(Op_Real_, "ComplexHypreParMatrix has no real part!");
+   return dynamic_cast<const HypreParMatrix &>(*Op_Real_);
+}
+
+const HypreParMatrix & ComplexHypreParMatrix::imag() const
+{
+   MFEM_ASSERT(Op_Imag_, "ComplexHypreParMatrix has no imaginary part!");
+   return dynamic_cast<const HypreParMatrix &>(*Op_Imag_);
+}
+
+HypreParMatrix * ComplexHypreParMatrix::GetSystemMatrix() const
+{
+   HypreParMatrix * A_r = dynamic_cast<HypreParMatrix*>(Op_Real_);
+   HypreParMatrix * A_i = dynamic_cast<HypreParMatrix*>(Op_Imag_);
+
+   if ( A_r == NULL && A_i == NULL ) { return NULL; }
+
+   HYPRE_Int global_num_rows_r = (A_r) ? A_r->GetGlobalNumRows() : 0;
+   HYPRE_Int global_num_rows_i = (A_i) ? A_i->GetGlobalNumRows() : 0;
+   HYPRE_Int global_num_rows = std::max(global_num_rows_r, global_num_rows_i);
+
+   HYPRE_Int global_num_cols_r = (A_r) ? A_r->GetGlobalNumCols() : 0;
+   HYPRE_Int global_num_cols_i = (A_i) ? A_i->GetGlobalNumCols() : 0;
+   HYPRE_Int global_num_cols = std::max(global_num_cols_r, global_num_cols_i);
+
+   int row_starts_size = (HYPRE_AssumedPartitionCheck()) ? 2 : nranks_ + 1;
+   HYPRE_Int * row_starts = hypre_CTAlloc(HYPRE_Int, row_starts_size);
+   HYPRE_Int * col_starts = hypre_CTAlloc(HYPRE_Int, row_starts_size);
+
+   const HYPRE_Int * row_starts_z = (A_r) ? A_r->RowPart() :
+                                    ((A_i) ? A_i->RowPart() : NULL);
+   const HYPRE_Int * col_starts_z = (A_r) ? A_r->ColPart() :
+                                    ((A_i) ? A_i->ColPart() : NULL);
+
+   for (int i = 0; i < row_starts_size; i++)
+   {
+      row_starts[i] = 2 * row_starts_z[i];
+      col_starts[i] = 2 * col_starts_z[i];
+   }
+
+   SparseMatrix diag_r, diag_i, offd_r, offd_i;
+   HYPRE_Int * cmap_r, * cmap_i;
+
+   int nrows_r = 0, nrows_i = 0, ncols_r = 0, ncols_i = 0;
+   int ncols_offd_r = 0, ncols_offd_i = 0;
+   if (A_r)
+   {
+      A_r->GetDiag(diag_r);
+      A_r->GetOffd(offd_r, cmap_r);
+      nrows_r = diag_r.Height();
+      ncols_r = diag_r.Width();
+      ncols_offd_r = offd_r.Width();
+   }
+   if (A_i)
+   {
+      A_i->GetDiag(diag_i);
+      A_i->GetOffd(offd_i, cmap_i);
+      nrows_i = diag_i.Height();
+      ncols_i = diag_i.Width();
+      ncols_offd_i = offd_i.Width();
+   }
+   int nrows = std::max(nrows_r, nrows_i);
+   int ncols = std::max(ncols_r, ncols_i);
+
+   // Determine the unique set of off-diagonal columns global indices
+   std::set<int> cset;
+   for (int i=0; i<ncols_offd_r; i++)
+   {
+      cset.insert(cmap_r[i]);
+   }
+   for (int i=0; i<ncols_offd_i; i++)
+   {
+      cset.insert(cmap_i[i]);
+   }
+   int num_cols_offd = (int)cset.size();
+
+   // Exatract pointers to the various CSR arrays of the diagonal blocks
+   const int * diag_r_I = (A_r) ? diag_r.GetI() : NULL;
+   const int * diag_i_I = (A_i) ? diag_i.GetI() : NULL;
+
+   const int * diag_r_J = (A_r) ? diag_r.GetJ() : NULL;
+   const int * diag_i_J = (A_i) ? diag_i.GetJ() : NULL;
+
+   const double * diag_r_D = (A_r) ? diag_r.GetData() : NULL;
+   const double * diag_i_D = (A_i) ? diag_i.GetData() : NULL;
+
+   int diag_r_nnz = (diag_r_I) ? diag_r_I[nrows] : 0;
+   int diag_i_nnz = (diag_i_I) ? diag_i_I[nrows] : 0;
+   int diag_nnz = 2 * (diag_r_nnz + diag_i_nnz);
+
+   // Exatract pointers to the various CSR arrays of the off-diagonal blocks
+   const int * offd_r_I = (A_r) ? offd_r.GetI() : NULL;
+   const int * offd_i_I = (A_i) ? offd_i.GetI() : NULL;
+
+   const int * offd_r_J = (A_r) ? offd_r.GetJ() : NULL;
+   const int * offd_i_J = (A_i) ? offd_i.GetJ() : NULL;
+
+   const double * offd_r_D = (A_r) ? offd_r.GetData() : NULL;
+   const double * offd_i_D = (A_i) ? offd_i.GetData() : NULL;
+
+   int offd_r_nnz = (offd_r_I) ? offd_r_I[nrows] : 0;
+   int offd_i_nnz = (offd_i_I) ? offd_i_I[nrows] : 0;
+   int offd_nnz = 2 * (offd_r_nnz + offd_i_nnz);
+
+   // Allocate CSR arrays for the combined matrix
+   HYPRE_Int * diag_I = hypre_CTAlloc(HYPRE_Int, 2 * nrows + 1);
+   HYPRE_Int * diag_J = hypre_CTAlloc(HYPRE_Int, diag_nnz);
+   double    * diag_D = hypre_CTAlloc(double, diag_nnz);
+
+   HYPRE_Int * offd_I = hypre_CTAlloc(HYPRE_Int, 2 * nrows + 1);
+   HYPRE_Int * offd_J = hypre_CTAlloc(HYPRE_Int, offd_nnz);
+   double    * offd_D = hypre_CTAlloc(double, offd_nnz);
+   HYPRE_Int * cmap   = hypre_CTAlloc(HYPRE_Int, 2 * num_cols_offd);
+
+   // Fill the CSR arrays for the diagonal portion of the matrix
+   const double factor = (convention_ == HERMITIAN) ? 1.0 : -1.0;
+
+   diag_I[0] = 0;
+   diag_I[nrows] = diag_r_nnz + diag_i_nnz;
+   for (int i=0; i<nrows; i++)
+   {
+      diag_I[i + 1]         = ((diag_r_I)?diag_r_I[i+1]:0) +
+                              ((diag_i_I)?diag_i_I[i+1]:0);
+      diag_I[i + nrows + 1] = diag_I[i+1] + diag_r_nnz + diag_i_nnz;
+
+      if (diag_r_I)
+      {
+         for (int j=0; j<diag_r_I[i+1] - diag_r_I[i]; j++)
+         {
+            diag_J[diag_I[i] + j] = diag_r_J[diag_r_I[i] + j];
+            diag_D[diag_I[i] + j] = diag_r_D[diag_r_I[i] + j];
+
+            diag_J[diag_I[i+nrows] + j] =
+               diag_r_J[diag_r_I[i] + j] + ncols;
+            diag_D[diag_I[i+nrows] + j] =
+               factor * diag_r_D[diag_r_I[i] + j];
+         }
+      }
+      if (diag_i_I)
+      {
+         const int off_r = (diag_r_I)?(diag_r_I[i+1] - diag_r_I[i]):0;
+         for (int j=0; j<diag_i_I[i+1] - diag_i_I[i]; j++)
+         {
+            diag_J[diag_I[i] + off_r + j] =  diag_i_J[diag_i_I[i] + j] + ncols;
+            diag_D[diag_I[i] + off_r + j] = -diag_i_D[diag_i_I[i] + j];
+
+            diag_J[diag_I[i+nrows] + off_r + j] = diag_i_J[diag_i_I[i] + j];
+            diag_D[diag_I[i+nrows] + off_r + j] =
+               factor * diag_i_D[diag_i_I[i] + j];
+         }
+      }
+   }
+
+   // Determine the mappings describing the layout of off-diagonal columns
+   int num_recv_procs = 0;
+   HYPRE_Int * offd_col_start_stop = NULL;
+   this->getColStartStop(A_r, A_i, num_recv_procs, offd_col_start_stop);
+
+   std::set<int>::iterator sit;
+   std::map<int,int> cmapa, cmapb, cinvmap;
+   for (sit=cset.begin(); sit!=cset.end(); sit++)
+   {
+      int col_orig = *sit;
+      int col_2x2  = -1;
+      int col_size = 0;
+      for (int i=0; i<num_recv_procs; i++)
+      {
+         if (offd_col_start_stop[2*i] <= col_orig &&
+             col_orig < offd_col_start_stop[2*i+1])
+         {
+            col_2x2 = offd_col_start_stop[2*i] + col_orig;
+            col_size = offd_col_start_stop[2*i+1] - offd_col_start_stop[2*i];
+            break;
+         }
+      }
+      cmapa[*sit] = col_2x2;
+      cmapb[*sit] = col_2x2 + col_size;
+      cinvmap[col_2x2] = -1;
+      cinvmap[col_2x2 + col_size] = -1;
+   }
+   delete [] offd_col_start_stop;
+
+   std::map<int, int>::iterator mit;
+   int i = 0;
+   for (mit=cinvmap.begin(); mit!=cinvmap.end(); mit++, i++)
+   {
+      mit->second = i;
+      cmap[i] = mit->first;
+   }
+
+   // Fill the CSR arrays for the off-diagonal portion of the matrix
+   offd_I[0] = 0;
+   offd_I[nrows] = offd_r_nnz + offd_i_nnz;
+   for (int i=0; i<nrows; i++)
+   {
+      offd_I[i + 1]         = ((offd_r_I)?offd_r_I[i+1]:0) +
+                              ((offd_i_I)?offd_i_I[i+1]:0);
+      offd_I[i + nrows + 1] = offd_I[i+1] + offd_r_nnz + offd_i_nnz;
+
+      if (offd_r_I)
+      {
+         const int off_i = (offd_i_I)?(offd_i_I[i+1] - offd_i_I[i]):0;
+         for (int j=0; j<offd_r_I[i+1] - offd_r_I[i]; j++)
+         {
+            offd_J[offd_I[i] + j] =
+               cinvmap[cmapa[cmap_r[offd_r_J[offd_r_I[i] + j]]]];
+            offd_D[offd_I[i] + j] = offd_r_D[offd_r_I[i] + j];
+
+            offd_J[offd_I[i+nrows] + off_i + j] =
+               cinvmap[cmapb[cmap_r[offd_r_J[offd_r_I[i] + j]]]];
+            offd_D[offd_I[i+nrows] + off_i + j] =
+               factor * offd_r_D[offd_r_I[i] + j];
+         }
+      }
+      if (offd_i_I)
+      {
+         const int off_r = (offd_r_I)?(offd_r_I[i+1] - offd_r_I[i]):0;
+         for (int j=0; j<offd_i_I[i+1] - offd_i_I[i]; j++)
+         {
+            offd_J[offd_I[i] + off_r + j] =
+               cinvmap[cmapb[cmap_i[offd_i_J[offd_i_I[i] + j]]]];
+            offd_D[offd_I[i] + off_r + j] = -offd_i_D[offd_i_I[i] + j];
+
+            offd_J[offd_I[i+nrows] + j] =
+               cinvmap[cmapa[cmap_i[offd_i_J[offd_i_I[i] + j]]]];
+            offd_D[offd_I[i+nrows] + j] = factor * offd_i_D[offd_i_I[i] + j];
+         }
+      }
+   }
+
+   // Construct the combined matrix
+   HypreParMatrix * A = new HypreParMatrix(comm_,
+                                           2 * global_num_rows,
+                                           2 * global_num_cols,
+                                           row_starts, col_starts,
+                                           diag_I, diag_J, diag_D,
+                                           offd_I, offd_J, offd_D,
+                                           2 * num_cols_offd, cmap);
+
+   // Give the new matrix ownership of its interanl arrays
+   A->SetOwnerFlags(-1,-1,-1);
+   hypre_CSRMatrixSetDataOwner(((hypre_ParCSRMatrix*)(*A))->diag,1);
+   hypre_CSRMatrixSetDataOwner(((hypre_ParCSRMatrix*)(*A))->offd,1);
+   hypre_ParCSRMatrixSetRowStartsOwner((hypre_ParCSRMatrix*)(*A),1);
+   hypre_ParCSRMatrixSetColStartsOwner((hypre_ParCSRMatrix*)(*A),1);
+
+   return A;
+}
+
+void
+ComplexHypreParMatrix::getColStartStop(const HypreParMatrix * A_r,
+                                       const HypreParMatrix * A_i,
+                                       int & num_recv_procs,
+                                       HYPRE_Int *& offd_col_start_stop) const
+{
+   hypre_ParCSRCommPkg * comm_pkg_r =
+      (A_r) ? hypre_ParCSRMatrixCommPkg((hypre_ParCSRMatrix*)(*A_r)) : NULL;
+   hypre_ParCSRCommPkg * comm_pkg_i =
+      (A_i) ? hypre_ParCSRMatrixCommPkg((hypre_ParCSRMatrix*)(*A_i)) : NULL;
+
+   std::set<HYPRE_Int> send_procs, recv_procs;
+   if ( comm_pkg_r )
+   {
+      for (HYPRE_Int i=0; i<comm_pkg_r->num_sends; i++)
+      {
+         send_procs.insert(comm_pkg_r->send_procs[i]);
+      }
+      for (HYPRE_Int i=0; i<comm_pkg_r->num_recvs; i++)
+      {
+         recv_procs.insert(comm_pkg_r->recv_procs[i]);
+      }
+   }
+   if ( comm_pkg_i )
+   {
+      for (HYPRE_Int i=0; i<comm_pkg_i->num_sends; i++)
+      {
+         send_procs.insert(comm_pkg_i->send_procs[i]);
+      }
+      for (HYPRE_Int i=0; i<comm_pkg_i->num_recvs; i++)
+      {
+         recv_procs.insert(comm_pkg_i->recv_procs[i]);
+      }
+   }
+
+   num_recv_procs = (int)recv_procs.size();
+
+   HYPRE_Int loc_start_stop[2];
+   offd_col_start_stop = new HYPRE_Int[2 * num_recv_procs];
+
+   const HYPRE_Int * row_part = (A_r) ? A_r->RowPart() :
+                                ((A_i) ? A_i->RowPart() : NULL);
+
+   int row_part_ind = (HYPRE_AssumedPartitionCheck()) ? 0 : myid_;
+   loc_start_stop[0] = row_part[row_part_ind];
+   loc_start_stop[1] = row_part[row_part_ind+1];
+
+   MPI_Request * req = new MPI_Request[send_procs.size()+recv_procs.size()];
+   MPI_Status * stat = new MPI_Status[send_procs.size()+recv_procs.size()];
+   int send_count = 0;
+   int recv_count = 0;
+   int tag = 0;
+
+   std::set<HYPRE_Int>::iterator sit;
+   for (sit=send_procs.begin(); sit!=send_procs.end(); sit++)
+   {
+      MPI_Isend(loc_start_stop, 2, HYPRE_MPI_INT,
+                *sit, tag, comm_, &req[send_count]);
+      send_count++;
+   }
+   for (sit=recv_procs.begin(); sit!=recv_procs.end(); sit++)
+   {
+      MPI_Irecv(&offd_col_start_stop[2*recv_count], 2, HYPRE_MPI_INT,
+                *sit, tag, comm_, &req[send_count+recv_count]);
+      recv_count++;
+   }
+
+   MPI_Waitall(send_count+recv_count, req, stat);
+
+   delete [] req;
+   delete [] stat;
+}
+
+
+#endif // MFEM_USE_MPI
+
 }

linalg/complex_operator.hpp

@@ -14,6 +14,9 @@
 
 #include "operator.hpp"
 #include "sparsemat.hpp"
+#ifdef MFEM_USE_MPI
+#include "hypre.hpp"
+#endif
 
 namespace mfem
 {
@@ -27,7 +30,8 @@ namespace mfem
     ComplexOperator allows one to choose a convention upon construction, which
     facilitates symmetry.
 
-    Matrix-vector products are then computed as:
+    If we let (y_r + i y_i) = (Op_r + i Op_i)(x_r + i x_i) then Matrix-vector
+    products are then computed as:
 
     1. When Convention::HERMITIAN is used (default)
     / y_r \   / Op_r -Op_i \ / x_r \
@@ -38,6 +42,8 @@ namespace mfem
     / y_r \   / Op_r -Op_i \ / x_r \
     |     | = |            | |     |
     \-y_i /   \-Op_i -Op_r / \ x_i /
+    In other words, Matrix-vector products with Convention::BLOCK_SYMMETRIC
+    compute the complex conjugate of Op*x.
 
     Either convention can be used with a given complex operator,
     however, each of them is best suited for certain classes of
@@ -82,9 +88,30 @@ public:
 
    virtual ~ComplexOperator();
 
+   /** @brief Check for existence of real or imaginary part of the operator
+
+       These methods do not check that the operators are non-zero but
+       only that the operators have been set.
+    */
+   bool hasRealPart() const { return Op_Real_ != NULL; }
+   bool hasImagPart() const { return Op_Imag_ != NULL; }
+
+   /** @brief Real or imaginary part accessor methods
+
+       The following accessor methods should only be called if the
+       requested part of the opertor is known to exist.  This
+       can be checked with hasRealPart() or hasImagPart().
+   */
+   virtual Operator & real();
+   virtual Operator & imag();
+   virtual const Operator & real() const;
+   virtual const Operator & imag() const;
+
    virtual void Mult(const Vector &x, Vector &y) const;
    virtual void MultTranspose(const Vector &x, Vector &y) const;
 
+   virtual Type GetType() const { return Complex_Operator; }
+
 protected:
    // Let this be hidden from the public interface since the implementation
    // depends on internal members
@@ -127,9 +154,74 @@ public:
       : ComplexOperator(A_Real, A_Imag, ownReal, ownImag, convention)
    {}
 
+   virtual SparseMatrix & real();
+   virtual SparseMatrix & imag();
+
+   virtual const SparseMatrix & real() const;
+   virtual const SparseMatrix & imag() const;
+
+   /** Combine the blocks making up this complex operator into a
+       single SparseMatrix.  The resulting matrix can be passed to
+       solvers which require access to the matrix entries themselves,
+       such as sparse direct solvers, rather than simply the action of
+       the opertor.  Note that this combined operator requires roughly
+       twice the memory of the block structured operator. */
    SparseMatrix * GetSystemMatrix() const;
+
+   virtual Type GetType() const { return MFEM_ComplexSparseMat; }
 };
 
+#ifdef MFEM_USE_MPI
+
+/** @brief Specialization of the ComplexOperator built from a pair of
+    HypreParMatrices.
+
+    The purpose of this specialization is to construct a single
+    HypreParMatrix object which is equivalent to the 2x2 block system
+    that the ComplexOperator mimics. The resulting HypreParMatrix can
+    then be passed along to solvers which require access to the CSR
+    matrix data such as SuperLU, STRUMPACK, or similar sparse linear
+    solvers.
+
+    See ComplexOperator documentation in operator.hpp for more information.
+ */
+class ComplexHypreParMatrix : public ComplexOperator
+{
+public:
+   ComplexHypreParMatrix(HypreParMatrix * A_Real, HypreParMatrix * A_Imag,
+                         bool ownReal, bool ownImag,
+                         Convention convention = HERMITIAN);
+
+   virtual HypreParMatrix & real();
+   virtual HypreParMatrix & imag();
+
+   virtual const HypreParMatrix & real() const;
+   virtual const HypreParMatrix & imag() const;
+
+   /** Combine the blocks making up this complex operator into a
+       single HypreParMatrix.  The resulting matrix can be passed to
+       solvers which require access to the matrix entries themselves,
+       such as sparse direct solvers or Hypre preconditioners, rather
+       than simply the action of the opertor.  Note that this combined
+       operator requires roughly twice the memory of the block
+       structured operator. */
+   HypreParMatrix * GetSystemMatrix() const;
+
+   virtual Type GetType() const { return Complex_Hypre_ParCSR; }
+
+private:
+   void getColStartStop(const HypreParMatrix * A_r,
+                        const HypreParMatrix * A_i,
+                        int & num_recv_procs,
+                        HYPRE_Int *& offd_col_start_stop) const;
+
+   MPI_Comm comm_;
+   int myid_;
+   int nranks_;
+};
+
+#endif // MFEM_USE_MPI
+
 }
 
-#endif
+#endif // MFEM_COMPLEX_OPERATOR

linalg/operator.hpp

@@ -124,14 +124,17 @@ public:
    enum Type
    {
       ANY_TYPE,         ///< ID for the base class Operator, i.e. any type.
-      MFEM_SPARSEMAT,   ///< ID for class SparseMatrix
+      MFEM_SPARSEMAT,   ///< ID for class SparseMatrix.
       Hypre_ParCSR,     ///< ID for class HypreParMatrix.
       PETSC_MATAIJ,     ///< ID for class PetscParMatrix, MATAIJ format.
       PETSC_MATIS,      ///< ID for class PetscParMatrix, MATIS format.
       PETSC_MATSHELL,   ///< ID for class PetscParMatrix, MATSHELL format.
       PETSC_MATNEST,    ///< ID for class PetscParMatrix, MATNEST format.
       PETSC_MATHYPRE,   ///< ID for class PetscParMatrix, MATHYPRE format.
-      PETSC_MATGENERIC  ///< ID for class PetscParMatrix, unspecified format.
+      PETSC_MATGENERIC, ///< ID for class PetscParMatrix, unspecified format.
+      Complex_Operator, ///< ID for class ComplexOperator.
+      MFEM_ComplexSparseMat, ///< ID for class ComplexSparseMatrix.
+      Complex_Hypre_ParCSR   ///< ID for class ComplexHypreParMatrix.
    };
 
    /// Return the type ID of the Operator class.
@@ -302,6 +305,24 @@ public:
 };
 
 
+/// Scaled Operator B: x -> a A(x).
+class ScaledOperator : public Operator
+{
+private:
+   const Operator &A_;
+   double a_;
+
+public:
+   /// Create a scalar product operator related to A.
+   explicit ScaledOperator(const Operator *A, double a)
+      : Operator(A->Width(), A->Height()), A_(*A), a_(a) { }
+
+   /// Operator application
+   virtual void Mult(const Vector &x, Vector &y) const
+   { A_.Mult(x, y); y *= a_; }
+};
+
+
 /** @brief The transpose of a given operator. Switches the roles of the methods
     Mult() and MultTranspose(). */
 class TransposeOperator : public Operator