Merge branch 'master' into gf-conv-dev

# Conflicts:
#	fem/fe.cpp

examples/ex_sconv_p.cpp

@@ -0,0 +1,503 @@
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+static double freq = 0.5, kappa;
+static int dim;
+
+double u_func(const Vector &);
+
+enum SCA_TYPE {INVALID_SCA_TYPE = -1,
+               H1_TYPE = 0,
+               L2_TYPE,
+               L2I_TYPE,
+               NUM_SCA_TYPES
+              };
+enum CONV_TYPE {INVALID_CONV_TYPE = -1,
+                PROJECTION = 0,
+                INTERPOLATION_OP,
+                SOLVE,
+                SOLVE_W_DBC,
+                NUM_CONV_TYPES
+               };
+
+FiniteElementCollection * GetFECollection(SCA_TYPE type, int p);
+ParFiniteElementSpace * GetFESpace(SCA_TYPE type, ParMesh &pmesh,
+                                   FiniteElementCollection &fec);
+string GetTypeName(SCA_TYPE type);
+string GetConvTypeName(CONV_TYPE type);
+string GetConvTypeShortName(CONV_TYPE type);
+
+void Projection(const ParGridFunction &v0, ParGridFunction &v1);
+void InterpolationOp(const ParGridFunction &v0, ParGridFunction &v1);
+void LeastSquares(SCA_TYPE t0, const ParGridFunction &v0,
+                  SCA_TYPE t1, ParGridFunction &v1);
+void LeastSquaresBC(SCA_TYPE t0, const ParGridFunction &v0,
+                    SCA_TYPE t1, ParGridFunction &v1,
+                    Coefficient &c);
+
+int main(int argc, char *argv[])
+{
+   // 1. Initialize MPI.
+   MPI_Session mpi(argc, argv);
+
+   // 2. Parse command-line options.
+   const char *mesh_file = "../data/star.mesh";
+   int ser_ref_levels = 0;
+   int par_ref_levels = 0;
+   int order0 = 1;
+   int order1 = 1;
+   int type0 = 0;
+   int type1 = 1;
+   int conv_type = -1;
+   bool static_cond = false;
+   bool pa = false;
+   const char *device_config = "cpu";
+   bool visualization = 1;
+
+   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(&order0, "-o0", "--initial-order",
+                  "Finite element order (polynomial degree) "
+                  "for initial field.");
+   args.AddOption(&order1, "-o1", "--final-order",
+                  "Finite element order (polynomial degree) "
+                  "for final field.");
+   args.AddOption(&type0, "-t0", "--initial-type",
+                  "Set the basis type for the initial field: "
+                  "0-H1, 1-L2, 2-L2I, -1 loop over all.");
+   args.AddOption(&type1, "-t1", "--final-type",
+                  "Set the basis type for the final field: "
+                  "0-H1, 1-L2, 2-L2I, -1 loop over all.");
+   args.AddOption(&conv_type, "-c", "--conversion-type",
+                  "Set the conversion scheme: "
+                  "0-Projection, 1-Interpolation Op, 2-Least Squares, "
+                  "3-Least Squares with BC, -1 loop over all.");
+   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(&pa, "-pa", "--partial-assembly", "-no-pa",
+                  "--no-partial-assembly", "Enable Partial Assembly.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+
+   args.Parse();
+   if (!args.Good())
+   {
+      if (mpi.Root()) { args.PrintUsage(cout); }
+      return 1;
+   }
+   if (mpi.Root()) { args.PrintOptions(cout); }
+   kappa = freq * M_PI;
+
+   // 3. Enable hardware devices such as GPUs, and programming models such as
+   //    CUDA, OCCA, RAJA and OpenMP based on command line options.
+   Device device(device_config);
+   if (mpi.Root()) { device.Print(); }
+
+   // 4. 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();
+
+   // 5. Refine the serial mesh on all processors to increase the resolution. In
+   //    this example we do 'ref_levels' of uniform refinement (2 by default, or
+   //    specified on the command line with -rs).
+   for (int lev = 0; lev < ser_ref_levels; lev++)
+   {
+      mesh->UniformRefinement();
+   }
+
+   // 6. Define a parallel mesh by a partitioning of the serial mesh. Refine
+   //    this mesh further in parallel to increase the resolution (1 time by
+   //    default, or specified on the command line with -rp). Once the parallel
+   //    mesh is defined, the serial mesh can be deleted.
+   ParMesh pmesh(MPI_COMM_WORLD, *mesh);
+   delete mesh;
+   for (int lev = 0; lev < par_ref_levels; lev++)
+   {
+      pmesh.UniformRefinement();
+   }
+
+   FunctionCoefficient uCoef(u_func);
+
+   int Ww = 300, Wh = 220, Fw = 3, Fh = 23, Ws = 15;
+
+   if (mpi.Root())
+   {
+      cout << "L2 Errors:" << endl;
+   }
+   int t0a = (type0 == -1) ? 0 : type0;
+   int t0b = (type0 == -1) ? NUM_SCA_TYPES : (type0+1);
+   for (int t0 = t0a; t0 < t0b; t0++)
+   {
+      FiniteElementCollection *fec0 = GetFECollection((SCA_TYPE)t0, order0);
+      ParFiniteElementSpace   *fes0 = GetFESpace((SCA_TYPE)t0, pmesh, *fec0);
+
+      ParGridFunction x0(fes0);
+      x0.ProjectCoefficient(uCoef);
+
+      double err0 = x0.ComputeL2Error(uCoef);
+      if (mpi.Root())
+      {
+         cout << "Initial " << GetTypeName((SCA_TYPE)t0)
+              << ": \t\t" << err0 << endl;
+      }
+
+      // nn. Send the solution by socket to a GLVis server.
+      if (visualization)
+      {
+         ostringstream oss;
+         oss << GetTypeName((SCA_TYPE)t0) << "(" << order0 << ")";
+         char vishost[] = "localhost";
+         int  visport   = 19916;
+         socketstream sol_sock0(vishost, visport);
+         sol_sock0 << "parallel " << pmesh.GetNRanks() << ' '
+                   << pmesh.GetMyRank() << '\n';
+         sol_sock0.precision(8);
+         sol_sock0 << "solution\n" << pmesh << x0
+                   << "window_title '" << oss.str() << "'"
+                   << "window_geometry "
+                   << Ws * (t0 - t0a) << " " << Ws * (t0 - t0a) << " "
+                   << (int)(1.5 * Ww) << " " << (int)(1.5 * Wh)
+                   << flush;
+      }
+
+      int t1a = (type1 == -1) ? 0 : type1;
+      int t1b = (type1 == -1) ? NUM_SCA_TYPES : (type1+1);
+      for (int t1 = t1a; t1 < t1b; t1++)
+      {
+         FiniteElementCollection *fec1 = GetFECollection((SCA_TYPE)t1, order1);
+         ParFiniteElementSpace   *fes1 = GetFESpace((SCA_TYPE)t1, pmesh, *fec1);
+
+         ParGridFunction y1(fes1);
+
+         if (mpi.Root())
+         {
+            cout << GetTypeName((SCA_TYPE)t0) << "(" << order0 << ")"
+                 << " -> "
+                 << GetTypeName((SCA_TYPE)t1) << "(" << order1 << ")"
+                 << ":" << endl;
+         }
+
+         int c01a = (conv_type == -1) ? 0 : conv_type;
+         int c01b = (conv_type == -1) ? NUM_CONV_TYPES : (conv_type+1);
+         for (int c01 = c01a; c01 < c01b; c01++)
+         {
+            string cmnt = "";
+
+            switch ((CONV_TYPE)c01)
+            {
+               case PROJECTION:
+                  Projection(x0, y1);
+                  break;
+               case INTERPOLATION_OP:
+                  cmnt = (t0 == (int)H1_TYPE) || (t0 == t1) ?
+                         "(should match projection)" : "(not expected to succeed)";
+                  InterpolationOp(x0, y1);
+                  break;
+               case SOLVE:
+                  LeastSquares((SCA_TYPE)t0, x0, (SCA_TYPE)t1, y1);
+                  break;
+               case SOLVE_W_DBC:
+                  LeastSquaresBC((SCA_TYPE)t0, x0, (SCA_TYPE)t1, y1, uCoef);
+                  break;
+               default:
+                  y1 = 0.0;
+            }
+
+            double err1 = y1.ComputeL2Error(uCoef);
+            cout << GetConvTypeName((CONV_TYPE)c01)
+                 << "\t\t" << err1 << "\t" << cmnt << endl;
+
+            if (visualization)
+            {
+               ostringstream oss;
+               oss << GetTypeName((SCA_TYPE)t0) << "(" << order0 << ")" << " --"
+                   << GetConvTypeShortName((CONV_TYPE)c01) << "--> "
+                   << GetTypeName((SCA_TYPE)t1)<< "(" << order1 << ")";
+               char vishost[] = "localhost";
+               int  visport   = 19916;
+               socketstream sol_sock1(vishost, visport);
+               sol_sock1 << "parallel " << pmesh.GetNRanks() << ' '
+                         << pmesh.GetMyRank() << '\n';
+               sol_sock1.precision(8);
+               sol_sock1 << "solution\n" << pmesh << y1
+                         << "window_title '" << oss.str() << "'"
+                         << "window_geometry "
+                         << (int)((Ww + Fw) * (1.5 + c01 - c01a) +
+                                  Ws * (t0 - t0a))
+                         << " " << (Wh + Fh) * (t1 - t1a) + Ws * (t0 - t0a)
+                         << " " << Ww << " " << Wh
+                         << flush;
+            }
+         }
+         if (mpi.Root())
+         {
+            cout << endl;
+         }
+
+         delete fes1;
+         delete fec1;
+      }
+
+      delete fes0;
+      delete fec0;
+
+      if (t0 < t0b - 1)
+      {
+         char c;
+         if (mpi.Root())
+         {
+            cout << "press (q)uit or (c)ontinue --> " << flush;
+            cin >> c;
+         }
+         MPI_Bcast(&c, 1, MPI_CHAR, 0, MPI_COMM_WORLD);
+
+         if (c != 'c')
+         {
+            break;
+         }
+      }
+      if (mpi.Root())
+      {
+         cout << endl;
+      }
+   }
+
+   return 0;
+}
+
+double u_func(const Vector &x)
+{
+   double kx = kappa * x[0];
+   double ky = kappa * x[1];
+   double kz = (dim == 3) ? (kappa * x[2]) : 0.0;
+
+   // Add the gradient of a scalar function
+   return cos(kx) * cos(ky) * cos(kz);
+}
+
+FiniteElementCollection * GetFECollection(SCA_TYPE type, int p)
+{
+   switch (type)
+   {
+      case H1_TYPE:
+         return new H1_FECollection(p, dim);
+      case L2_TYPE:
+         return new L2_FECollection(p-1, dim);
+      case L2I_TYPE:
+         return new L2_FECollection(p-1, dim, BasisType::GaussLegendre,
+                                    FiniteElement::INTEGRAL);
+      default:
+         return NULL;
+   }
+}
+
+ParFiniteElementSpace * GetFESpace(SCA_TYPE type,
+                                   ParMesh &pmesh,
+                                   FiniteElementCollection &fec)
+{
+   return new ParFiniteElementSpace(&pmesh, &fec);
+}
+
+string GetTypeName(SCA_TYPE type)
+{
+   switch (type)
+   {
+      case H1_TYPE:
+         return "    H1";
+      case L2_TYPE:
+         return "    L2";
+      case L2I_TYPE:
+         return "    L2I";
+      default:
+         return "--";
+   }
+}
+
+string GetConvTypeName(CONV_TYPE type)
+{
+   switch (type)
+   {
+      case PROJECTION:
+         return "Projection            ";
+      case INTERPOLATION_OP:
+         return "Interpolation Operator";
+      case SOLVE:
+         return "Least Squares         ";
+      case SOLVE_W_DBC:
+         return "Least Squares with BC ";
+      default:
+         return "--";
+   }
+}
+
+string GetConvTypeShortName(CONV_TYPE type)
+{
+   switch (type)
+   {
+      case PROJECTION:
+         return "Proj";
+      case INTERPOLATION_OP:
+         return "Interp";
+      case SOLVE:
+         return "LS";
+      case SOLVE_W_DBC:
+         return "LSwBC";
+      default:
+         return "--";
+   }
+}
+
+/** Perform a naive projection from one scalar field to another.
+
+   This scheme simply evaluates v0 at the interpolation points of v1.
+
+   If v0 has reduced continuity compared to v1 this can produce
+   results that depend on the order in which the elements are
+   traversed.
+
+   Suitable conversions:
+   H1 -> L2
+   H1 -> DG (same as L2)
+ */
+void Projection(const ParGridFunction &v0, ParGridFunction &v1)
+{
+   GridFunctionCoefficient v0Coef(&v0);
+   v1.ProjectCoefficient(v0Coef);
+}
+
+/** In theory this interpolation scheme should be equivalent to projection.
+
+    Building an interpolastion matrix could lead to computational
+    efficiency compared to simple projection if the operator will be
+    used several times.
+
+    Unfortunately this is broken for several combinations of source
+    and target fields.
+*/
+void InterpolationOp(const ParGridFunction &v0, ParGridFunction &v1)
+{
+   ParDiscreteLinearOperator op(v0.ParFESpace(), v1.ParFESpace());
+   op.AddDomainInterpolator(new IdentityInterpolator);
+   op.Assemble();
+   op.Finalize();
+
+   op.Mult(v0, v1);
+}
+
+/** Compute a least-squares best fit using the target basis functions.
+
+    This scheme is more difficult to setup and more computationally
+    expensive but the results can be significantly better than simple
+    projections.
+*/
+void LeastSquares(SCA_TYPE t0, const ParGridFunction &v0,
+                  SCA_TYPE t1, ParGridFunction &v1)
+{
+   ParFiniteElementSpace *fes0, *fes1;
+   fes0 = v0.ParFESpace();
+   fes1 = v1.ParFESpace();
+
+   ParMixedBilinearForm op(fes0, fes1);
+   op.AddDomainIntegrator(new MassIntegrator);
+   op.Assemble();
+   op.Finalize();
+
+   ParLinearForm b(v1.ParFESpace());
+   op.Mult(v0, b);
+
+   ParBilinearForm m(v1.ParFESpace());
+   m.AddDomainIntegrator(new MassIntegrator);
+   m.Assemble();
+   m.Finalize();
+
+   HypreParMatrix * M = m.ParallelAssemble();
+
+   HypreDiagScale diag(*M);
+   HyprePCG pcg(*M);
+   pcg.SetPreconditioner(diag);
+   pcg.SetTol(1e-12);
+   pcg.SetMaxIter(1000);
+
+   Vector B, X;
+   b.ParallelAssemble(B);
+
+   X.SetSize(v1.ParFESpace()->TrueVSize()); X = 0.0;
+   pcg.Mult(B, X);
+   v1.Distribute(X);
+
+   delete M;
+}
+
+/** Compute a least-squares best fit with boundary conditions.
+
+    This scheme is virtually identical to the previous one but it
+    makes use of boundary values, when available, to improve the
+    accuracy.  This scheme can produce significantly better results
+    when the normal derivative of the field is large near the
+    boundary.  This is particularly true when the field is
+    under-resolved near the boundary.
+*/
+void LeastSquaresBC(SCA_TYPE t0, const ParGridFunction &v0,
+                    SCA_TYPE t1, ParGridFunction &v1,
+                    Coefficient &c)
+{
+   ParFiniteElementSpace *fes0, *fes1;
+   fes0 = v0.ParFESpace();
+   fes1 = v1.ParFESpace();
+
+   ParMixedBilinearForm op(fes0, fes1);
+   op.AddDomainIntegrator(new MassIntegrator);
+   op.Assemble();
+   op.Finalize();
+
+   ParLinearForm b(v1.ParFESpace());
+   op.Mult(v0, b);
+
+   ParBilinearForm m(v1.ParFESpace());
+   m.AddDomainIntegrator(new MassIntegrator);
+   m.Assemble();
+   m.Finalize();
+
+   Array<int> ess_bdr;
+   Array<int> ess_tdof_list;
+   if (v1.ParFESpace()->GetParMesh()->bdr_attributes.Size())
+   {
+      ess_bdr.SetSize(v1.ParFESpace()->GetParMesh()->bdr_attributes.Max());
+      ess_bdr = 1;
+      v1.ParFESpace()->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+
+   if (t1 == H1_TYPE)
+   {
+      v1.ProjectBdrCoefficient(c, ess_bdr);
+   }
+
+   OperatorPtr M;
+   Vector B, X;
+   m.FormLinearSystem(ess_tdof_list, v1, b, M, X, B);
+
+   HypreDiagScale diag(*M.As<HypreParMatrix>());
+   HyprePCG pcg(*M.As<HypreParMatrix>());
+   pcg.SetPreconditioner(diag);
+   pcg.SetTol(1e-12);
+   pcg.SetMaxIter(1000);
+
+   pcg.Mult(B, X);
+   v1.Distribute(X);
+}

examples/ex_vconv_p.cpp

@@ -0,0 +1,586 @@
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+static double freq = 0.5, kappa;
+static int dim;
+
+void u_func(const Vector &, Vector &);
+
+enum VEC_TYPE {INVALID_VEC_TYPE = -1,
+               H1V_TYPE = 0,
+               ND_TYPE,
+               RT_TYPE,
+               L2V_TYPE,
+               NUM_VEC_TYPES
+              };
+enum CONV_TYPE {INVALID_CONV_TYPE = -1,
+                PROJECTION = 0,
+                INTERPOLATION_OP,
+                SOLVE,
+                SOLVE_W_DBC,
+                NUM_CONV_TYPES
+               };
+
+FiniteElementCollection * GetFECollection(VEC_TYPE type, int p);
+ParFiniteElementSpace * GetFESpace(VEC_TYPE type, ParMesh &pmesh,
+                                   FiniteElementCollection &fec);
+string GetTypeName(VEC_TYPE type);
+string GetConvTypeName(CONV_TYPE type);
+
+void Projection(const ParGridFunction &v0, ParGridFunction &v1);
+void InterpolationOp(const ParGridFunction &v0, ParGridFunction &v1);
+void LeastSquares(VEC_TYPE t0, const ParGridFunction &v0,
+                  VEC_TYPE t1, ParGridFunction &v1);
+void LeastSquaresBC(VEC_TYPE t0, const ParGridFunction &v0,
+                    VEC_TYPE t1, ParGridFunction &v1,
+                    VectorCoefficient &vc);
+
+int main(int argc, char *argv[])
+{
+   // 1. Initialize MPI.
+   MPI_Session mpi(argc, argv);
+
+   // 2. Parse command-line options.
+   const char *mesh_file = "../data/star.mesh";
+   int ser_ref_levels = 0;
+   int par_ref_levels = 0;
+   int order0 = 1;
+   int order1 = 1;
+   int type0 = 0;
+   int type1 = 1;
+   int conv_type = -1;
+   bool static_cond = false;
+   bool pa = false;
+   const char *device_config = "cpu";
+   bool visualization = 1;
+
+   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(&order0, "-o0", "--initial-order",
+                  "Finite element order (polynomial degree) "
+                  "for initial field.");
+   args.AddOption(&order1, "-o1", "--final-order",
+                  "Finite element order (polynomial degree) "
+                  "for final field.");
+   args.AddOption(&type0, "-t0", "--initial-type",
+                  "Set the basis type for the initial field: "
+                  "0-H1V, 1-H(Curl), 2-H(Div), 3-L2V, -1 loop over all.");
+   args.AddOption(&type1, "-t1", "--final-type",
+                  "Set the basis type for the final field: "
+                  "0-H1V, 1-H(Curl), 2-H(Div), 3-L2V, -1 loop over all.");
+   args.AddOption(&conv_type, "-c", "--conversion-type",
+                  "Set the conversion scheme: "
+                  "0-Projection, 1-Interpolation Op, 2-Least Squares, "
+                  "3-Least Squares with BC, -1 loop over all.");
+   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(&pa, "-pa", "--partial-assembly", "-no-pa",
+                  "--no-partial-assembly", "Enable Partial Assembly.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+
+   args.Parse();
+   if (!args.Good())
+   {
+      if (mpi.Root()) { args.PrintUsage(cout); }
+      return 1;
+   }
+   if (mpi.Root()) { args.PrintOptions(cout); }
+   kappa = freq * M_PI;
+
+   // 3. Enable hardware devices such as GPUs, and programming models such as
+   //    CUDA, OCCA, RAJA and OpenMP based on command line options.
+   Device device(device_config);
+   if (mpi.Root()) { device.Print(); }
+
+   // 4. 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();
+
+   // 5. Refine the serial mesh on all processors to increase the resolution. In
+   //    this example we do 'ref_levels' of uniform refinement (2 by default, or
+   //    specified on the command line with -rs).
+   for (int lev = 0; lev < ser_ref_levels; lev++)
+   {
+      mesh->UniformRefinement();
+   }
+
+   // 6. Define a parallel mesh by a partitioning of the serial mesh. Refine
+   //    this mesh further in parallel to increase the resolution (1 time by
+   //    default, or specified on the command line with -rp). Once the parallel
+   //    mesh is defined, the serial mesh can be deleted.
+   ParMesh pmesh(MPI_COMM_WORLD, *mesh);
+   delete mesh;
+   for (int lev = 0; lev < par_ref_levels; lev++)
+   {
+      pmesh.UniformRefinement();
+   }
+
+   VectorFunctionCoefficient uCoef(dim, u_func);
+
+   int Ww = 300, Wh = 220, Fw = 3, Fh = 23, Ws = 15;
+
+   if (mpi.Root())
+   {
+      cout << "L2 Errors:" << endl;
+   }
+   int t0a = (type0 == -1) ? 0 : type0;
+   int t0b = (type0 == -1) ? NUM_VEC_TYPES : (type0+1);
+   for (int t0 = t0a; t0 < t0b; t0++)
+   {
+      FiniteElementCollection *fec0 = GetFECollection((VEC_TYPE)t0, order0);
+      ParFiniteElementSpace   *fes0 = GetFESpace((VEC_TYPE)t0, pmesh, *fec0);
+
+      ParGridFunction x0(fes0);
+      x0.ProjectCoefficient(uCoef);
+
+      double err0 = x0.ComputeL2Error(uCoef);
+      if (mpi.Root())
+      {
+         cout << "Initial " << GetTypeName((VEC_TYPE)t0)
+              << ": \t\t" << err0 << endl;
+      }
+
+      // nn. Send the solution by socket to a GLVis server.
+      if (visualization)
+      {
+         ostringstream oss;
+         oss << GetTypeName((VEC_TYPE)t0);
+         char vishost[] = "localhost";
+         int  visport   = 19916;
+         socketstream sol_sock0(vishost, visport);
+         sol_sock0 << "parallel " << pmesh.GetNRanks() << ' '
+                   << pmesh.GetMyRank() << '\n';
+         sol_sock0.precision(8);
+         sol_sock0 << "solution\n" << pmesh << x0
+                   << "keys vvv "
+                   << "window_title '" << oss.str() << "'"
+                   << "window_geometry "
+                   << Ws * (t0 - t0a) << " " << Ws * (t0 - t0a) << " "
+                   << (int)(1.5 * Ww) << " " << (int)(1.5 * Wh)
+                   << flush;
+      }
+
+      int t1a = (type1 == -1) ? 0 : type1;
+      int t1b = (type1 == -1) ? NUM_VEC_TYPES : (type1+1);
+      for (int t1 = t1a; t1 < t1b; t1++)
+      {
+         FiniteElementCollection *fec1 = GetFECollection((VEC_TYPE)t1, order1);
+         ParFiniteElementSpace   *fes1 = GetFESpace((VEC_TYPE)t1, pmesh, *fec1);
+
+         ParGridFunction x1(fes1);
+
+         if (mpi.Root())
+         {
+            cout << GetTypeName((VEC_TYPE)t0) << " -> "
+                 << GetTypeName((VEC_TYPE)t1) << ":" << endl;
+         }
+
+         int c01a = (conv_type == -1) ? 0 : conv_type;
+         int c01b = (conv_type == -1) ? NUM_CONV_TYPES : (conv_type+1);
+         for (int c01 = c01a; c01 < c01b; c01++)
+         {
+            switch ((CONV_TYPE)c01)
+            {
+               case PROJECTION:
+                  Projection(x0, x1);
+                  break;
+               case INTERPOLATION_OP:
+                  // InterpolationOp(x0, x1);
+                  x1 = 0.0;
+                  break;
+               case SOLVE:
+                  LeastSquares((VEC_TYPE)t0, x0, (VEC_TYPE)t1, x1);
+                  break;
+               case SOLVE_W_DBC:
+                  LeastSquaresBC((VEC_TYPE)t0, x0, (VEC_TYPE)t1, x1, uCoef);
+                  break;
+               default:
+                  x1 = 0.0;
+            }
+
+            double err1 = x1.ComputeL2Error(uCoef);
+            cout << GetConvTypeName((CONV_TYPE)c01)
+                 << "\t\t" << err1 << endl;
+
+            if (visualization)
+            {
+               ostringstream oss;
+               oss << GetTypeName((VEC_TYPE)t0) << " --" << c01 << "--> "
+                   << GetTypeName((VEC_TYPE)t1);
+               char vishost[] = "localhost";
+               int  visport   = 19916;
+               socketstream sol_sock1(vishost, visport);
+               sol_sock1 << "parallel " << pmesh.GetNRanks() << ' '
+                         << pmesh.GetMyRank() << '\n';
+               sol_sock1.precision(8);
+               sol_sock1 << "solution\n" << pmesh << x1
+                         << "keys vvv "
+                         << "window_title '" << oss.str() << "'"
+                         << "window_geometry "
+                         << (int)((Ww + Fw) * (1.5 + c01 - c01a) +
+                                  Ws * (t0 - t0a))
+                         << " " << (Wh + Fh) * (t1 - t1a) + Ws * (t0 - t0a)
+                         << " " << Ww << " " << Wh
+                         << flush;
+            }
+         }
+         if (mpi.Root())
+         {
+            cout << endl;
+         }
+
+         delete fes1;
+         delete fec1;
+      }
+
+      delete fes0;
+      delete fec0;
+
+      if (t0 < t0b - 1)
+      {
+         char c;
+         if (mpi.Root())
+         {
+            cout << "press (q)uit or (c)ontinue --> " << flush;
+            cin >> c;
+         }
+         MPI_Bcast(&c, 1, MPI_CHAR, 0, MPI_COMM_WORLD);
+
+         if (c != 'c')
+         {
+            break;
+         }
+      }
+      if (mpi.Root())
+      {
+         cout << endl;
+      }
+   }
+
+   return 0;
+}
+
+void u_func(const Vector &x, Vector &u)
+{
+   u.SetSize(dim);
+
+   double kx = kappa * x[0];
+   double ky = kappa * x[1];
+   double kz = (dim == 3) ? (kappa * x[2]) : 0.0;
+
+   // Add the gradient of a scalar function
+   u(0) = sin(kx) * cos(ky);
+   u(1) = cos(kx) * sin(ky);
+   if (dim == 3)
+   {
+      u(0) *= cos(kz);
+      u(1) *= cos(kz);
+      u(2) = cos(kx) * cos(ky) * sin(kz);
+   }
+
+   // Add the curl of a vector function
+   u(0) -= cos(kx) * sin(ky);
+   u(1) += sin(kx) * cos(ky);
+   if (dim == 3)
+   {
+      u(0) += cos(kx) * sin(kz);
+      u(1) -= cos(ky) * sin(kz);
+      u(2) += (sin(ky) - sin(kx)) * cos(kz);
+   }
+}
+
+FiniteElementCollection * GetFECollection(VEC_TYPE type, int p)
+{
+   switch (type)
+   {
+      case H1V_TYPE:
+         return new H1_FECollection(p, dim);
+      case ND_TYPE:
+         return new ND_FECollection(p, dim);
+      case RT_TYPE:
+         return new RT_FECollection(p-1, dim);
+      case L2V_TYPE:
+         return new L2_FECollection(p-1, dim);
+      default:
+         return NULL;
+   }
+}
+
+ParFiniteElementSpace * GetFESpace(VEC_TYPE type,
+                                   ParMesh &pmesh,
+                                   FiniteElementCollection &fec)
+{
+   switch (type)
+   {
+      case H1V_TYPE:
+      case L2V_TYPE:
+         return new ParFiniteElementSpace(&pmesh, &fec, dim);
+      case ND_TYPE:
+      case RT_TYPE:
+         return new ParFiniteElementSpace(&pmesh, &fec);
+      default:
+         return NULL;
+   }
+}
+
+string GetTypeName(VEC_TYPE type)
+{
+   switch (type)
+   {
+      case H1V_TYPE:
+         return "    H1V";
+      case ND_TYPE:
+         return "H(Curl)";
+      case RT_TYPE:
+         return " H(Div)";
+      case L2V_TYPE:
+         return "    L2V";
+      default:
+         return "--";
+   }
+}
+
+string GetConvTypeName(CONV_TYPE type)
+{
+   switch (type)
+   {
+      case PROJECTION:
+         return "Projection            ";
+      case INTERPOLATION_OP:
+         return "Interpolation Operator";
+      case SOLVE:
+         return "Least Squares         ";
+      case SOLVE_W_DBC:
+         return "Least Squares with BC ";
+      default:
+         return "--";
+   }
+}
+
+/** Perform a naive projection from one vector field to another.
+
+   This scheme simply evaluates v0 at the interpolation points of v1.
+
+   If v0 has reduced continuity compared to v1 this can produce
+   results that depend on the order in which the elements are
+   traversed.
+
+   Suitable conversions:
+   H1V     -> H(Curl), H(Div), or L2V
+   H(Curl) -> L2V
+   H(Div)  -> L2V
+
+ */
+void Projection(const ParGridFunction &v0, ParGridFunction &v1)
+{
+   VectorGridFunctionCoefficient v0Coef(&v0);
+   v1.ProjectCoefficient(v0Coef);
+}
+
+/** In theory this interpolation scheme should be equivalent to projection.
+
+    Building an interpolastion matrix could lead to computational
+    efficiency compared to simple projection if the operator will be
+    used several times.
+
+    Unfortunately this is broken for several combinations of source
+    and target fields.
+*/
+void InterpolationOp(const ParGridFunction &v0, ParGridFunction &v1)
+{
+   ParDiscreteLinearOperator op(v0.ParFESpace(), v1.ParFESpace());
+   op.AddDomainInterpolator(new IdentityInterpolator);
+   op.Assemble();
+   op.Finalize();
+
+   op.Mult(v0, v1);
+}
+
+/** Compute a least-squares best fit using the target basis functions.
+
+    This scheme is more difficult to setup and more computationally
+    expensive but the results can be significantly better than simple
+    projections.
+*/
+void LeastSquares(VEC_TYPE t0, const ParGridFunction &v0,
+                  VEC_TYPE t1, ParGridFunction &v1)
+{
+   bool trans = false;
+   ParFiniteElementSpace *fes0, *fes1;
+   if ((t0 == H1V_TYPE || t0 == L2V_TYPE) &&
+       (t1 == ND_TYPE || t1 == RT_TYPE))
+   {
+      fes0 = v1.ParFESpace();
+      fes1 = v0.ParFESpace();
+      trans = true;
+   }
+   else
+   {
+      fes0 = v0.ParFESpace();
+      fes1 = v1.ParFESpace();
+   }
+   ParMixedBilinearForm op(fes0, fes1);
+   if (t0 == ND_TYPE || t0 == RT_TYPE || t1 == ND_TYPE || t1 == RT_TYPE)
+   {
+      op.AddDomainIntegrator(new VectorFEMassIntegrator);
+   }
+   else
+   {
+      op.AddDomainIntegrator(new VectorMassIntegrator);
+   }
+   op.Assemble();
+   op.Finalize();
+
+   ParLinearForm b(v1.ParFESpace());
+   if (trans)
+   {
+      op.MultTranspose(v0, b);
+   }
+   else
+   {
+      op.Mult(v0, b);
+   }
+
+   ParBilinearForm m(v1.ParFESpace());
+   if (t1 == ND_TYPE || t1 == RT_TYPE)
+   {
+      m.AddDomainIntegrator(new VectorFEMassIntegrator);
+   }
+   else
+   {
+      m.AddDomainIntegrator(new VectorMassIntegrator);
+   }
+   m.Assemble();
+   m.Finalize();
+
+   HypreParMatrix * M = m.ParallelAssemble();
+
+   HypreDiagScale diag(*M);
+   HyprePCG pcg(*M);
+   pcg.SetPreconditioner(diag);
+   pcg.SetTol(1e-12);
+   pcg.SetMaxIter(1000);
+
+   Vector B, X;
+   b.ParallelAssemble(B);
+
+   X.SetSize(v1.ParFESpace()->TrueVSize()); X = 0.0;
+   pcg.Mult(B, X);
+   v1.Distribute(X);
+
+   delete M;
+}
+
+/** Compute a least-squares best fit with boundary conditions.
+
+    This scheme is virtually identical to the previous one but it
+    makes use of boundary values, when available, to improve the
+    accuracy.  This scheme can produce significantly better results
+    when the normal derivative of the field is large near the
+    boundary.  This is particularly true when the field is
+    under-resolved near the boundary.
+*/
+void LeastSquaresBC(VEC_TYPE t0, const ParGridFunction &v0,
+                    VEC_TYPE t1, ParGridFunction &v1,
+                    VectorCoefficient &vc)
+{
+   bool trans = false;
+   ParFiniteElementSpace *fes0, *fes1;
+   if ((t0 == H1V_TYPE || t0 == L2V_TYPE) &&
+       (t1 == ND_TYPE || t1 == RT_TYPE))
+   {
+      fes0 = v1.ParFESpace();
+      fes1 = v0.ParFESpace();
+      trans = true;
+   }
+   else
+   {
+      fes0 = v0.ParFESpace();
+      fes1 = v1.ParFESpace();
+   }
+   ParMixedBilinearForm op(fes0, fes1);
+   if (t0 == ND_TYPE || t0 == RT_TYPE || t1 == ND_TYPE || t1 == RT_TYPE)
+   {
+      op.AddDomainIntegrator(new VectorFEMassIntegrator);
+   }
+   else
+   {
+      op.AddDomainIntegrator(new VectorMassIntegrator);
+   }
+   op.Assemble();
+   op.Finalize();
+
+   ParLinearForm b(v1.ParFESpace());
+   if (trans)
+   {
+      op.MultTranspose(v0, b);
+   }
+   else
+   {
+      op.Mult(v0, b);
+   }
+
+   ParBilinearForm m(v1.ParFESpace());
+   if (t1 == ND_TYPE || t1 == RT_TYPE)
+   {
+      m.AddDomainIntegrator(new VectorFEMassIntegrator);
+   }
+   else
+   {
+      m.AddDomainIntegrator(new VectorMassIntegrator);
+   }
+   m.Assemble();
+   m.Finalize();
+
+   Array<int> ess_bdr;
+   Array<int> ess_tdof_list;
+   if (v1.ParFESpace()->GetParMesh()->bdr_attributes.Size())
+   {
+      ess_bdr.SetSize(v1.ParFESpace()->GetParMesh()->bdr_attributes.Max());
+      ess_bdr = 1;
+      v1.ParFESpace()->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+
+   if (t1 == H1V_TYPE)
+   {
+      v1.ProjectBdrCoefficient(vc, ess_bdr);
+   }
+   if (t1 == ND_TYPE)
+   {
+      v1.ProjectBdrCoefficientTangent(vc, ess_bdr);
+   }
+   else if (t1 == RT_TYPE)
+   {
+      v1.ProjectBdrCoefficientNormal(vc, ess_bdr);
+   }
+
+   OperatorPtr M;
+   Vector B, X;
+   m.FormLinearSystem(ess_tdof_list, v1, b, M, X, B);
+
+   HypreDiagScale diag(*M.As<HypreParMatrix>());
+   HyprePCG pcg(*M.As<HypreParMatrix>());
+   pcg.SetPreconditioner(diag);
+   pcg.SetTol(1e-12);
+   pcg.SetMaxIter(1000);
+
+   pcg.Mult(B, X);
+   v1.Distribute(X);
+}

general/array.cpp

@@ -0,0 +1,214 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+
+// Abstract array data type
+
+#include "array.hpp"
+#include "../general/forall.hpp"
+#include <fstream>
+#include <type_traits>
+
+namespace mfem
+{
+
+template <class T>
+void Array<T>::Print(std::ostream &os, int width) const
+{
+   for (int i = 0; i < size; i++)
+   {
+      os << data[i];
+      if ( !((i+1) % width) || i+1 == size )
+      {
+         os << '\n';
+      }
+      else
+      {
+         os << " ";
+      }
+   }
+}
+
+template <class T>
+void Array<T>::Save(std::ostream &os, int fmt) const
+{
+   if (fmt == 0)
+   {
+      os << size << '\n';
+   }
+   for (int i = 0; i < size; i++)
+   {
+      os << operator[](i) << '\n';
+   }
+}
+
+template <class T>
+void Array<T>::Load(std::istream &in, int fmt)
+{
+   if (fmt == 0)
+   {
+      int new_size;
+      in >> new_size;
+      SetSize(new_size);
+   }
+   for (int i = 0; i < size; i++)
+   {
+      in >> operator[](i);
+   }
+}
+
+template <class T>
+T Array<T>::Max() const
+{
+   MFEM_ASSERT(size > 0, "Array is empty with size " << size);
+
+   T max = operator[](0);
+   for (int i = 1; i < size; i++)
+   {
+      if (max < operator[](i))
+      {
+         max = operator[](i);
+      }
+   }
+
+   return max;
+}
+
+template <class T>
+T Array<T>::Min() const
+{
+   MFEM_ASSERT(size > 0, "Array is empty with size " << size);
+
+   T min = operator[](0);
+   for (int i = 1; i < size; i++)
+   {
+      if (operator[](i) < min)
+      {
+         min = operator[](i);
+      }
+   }
+
+   return min;
+}
+
+// Partial Sum
+template <class T>
+void Array<T>::PartialSum()
+{
+   T sum = static_cast<T>(0);
+   for (int i = 0; i < size; i++)
+   {
+      sum+=operator[](i);
+      operator[](i) = sum;
+   }
+}
+
+template <class T>
+void Array<T>::Abs()
+{
+   static_assert(std::is_arithmetic<T>::value, "Use with arithmetic types!");
+   const bool useDevice = UseDevice();
+   const int N = size;
+   auto y = ReadWrite(useDevice);
+   mfem::forall_switch(useDevice, N, [=] MFEM_HOST_DEVICE (int i)
+   {
+      y[i] = std::abs(y[i]);
+   });
+}
+
+// Sum
+template <class T>
+T Array<T>::Sum() const
+{
+   T sum = static_cast<T>(0);
+   for (int i = 0; i < size; i++)
+   {
+      sum+=operator[](i);
+   }
+
+   return sum;
+}
+
+template <class T>
+int Array<T>::IsSorted() const
+{
+   T val_prev = operator[](0), val;
+   for (int i = 1; i < size; i++)
+   {
+      val=operator[](i);
+      if (val < val_prev)
+      {
+         return 0;
+      }
+      val_prev = val;
+   }
+
+   return 1;
+}
+
+template <class T>
+bool Array<T>::IsConstant() const
+{
+   if (size < 2) { return true; }
+   const T v0 = data[0];
+   for (int i = 1; i < size; i++)
+   {
+      if (data[i] != v0)
+      {
+         return false;
+      }
+   }
+
+   return true;
+}
+
+template <class T>
+void Array2D<T>::Load(const char *filename, int fmt)
+{
+   std::ifstream in;
+   in.open(filename, std::ifstream::in);
+   MFEM_VERIFY(in.is_open(), "File " << filename << " does not exist.");
+   Load(in, fmt);
+   in.close();
+}
+
+template <class T>
+void Array2D<T>::Print(std::ostream &os, int width_)
+{
+   int height = this->NumRows();
+   int width  = this->NumCols();
+
+   for (int i = 0; i < height; i++)
+   {
+      os << "[row " << i << "]\n";
+      for (int j = 0; j < width; j++)
+      {
+         os << (*this)(i,j);
+         if ( (j+1) == width_ || (j+1) % width_ == 0 )
+         {
+            os << '\n';
+         }
+         else
+         {
+            os << ' ';
+         }
+      }
+   }
+}
+
+template class Array<char>;
+template class Array<int>;
+template class Array<long long>;
+template class Array<real_t>;
+template class Array2D<int>;
+template class Array2D<real_t>;
+
+} // namespace mfem

general/array.hpp

@@ -16,13 +16,9 @@
 #include "mem_manager.hpp"
 #include "device.hpp"
 #include "error.hpp"
-#include "forall.hpp"
 #include "globals.hpp"
-#include "reducers.hpp"
-#include "scan.hpp"
 
 #include <iostream>
-#include <fstream>
 #include <cstdlib>
 #include <cstring>
 #include <algorithm>
@@ -139,8 +135,6 @@ public:
    /// Return the device flag of the Memory object used by the Array
    bool UseDevice() const { return data.UseDevice(); }
 
-   void UseDevice(bool use_dev) { data.UseDevice(use_dev); }
-
    /// Return true if the data will be deleted by the Array
    inline bool OwnsData() const { return data.OwnsHostPtr(); }
 
@@ -281,11 +275,11 @@ public:
 
    /** @brief Find the maximal element in the array, using the comparison
        operator `<` for class T. */
-   inline T Max() const;
+   T Max() const;
 
    /** @brief Find the minimal element in the array, using the comparison
        operator `<` for class T. */
-   inline T Min() const;
+   T Min() const;
 
    /// Sorts the array in ascending order. This requires operator< to be defined for T.
    void Sort() { std::sort((T*)data, data + size); }
@@ -303,22 +297,22 @@ public:
    }
 
    /// Return 1 if the array is sorted from lowest to highest.  Otherwise return 0.
-   inline int IsSorted() const;
+   int IsSorted() const;
 
    /// Does the Array have Size zero.
    bool IsEmpty() const { return Size() == 0; }
 
    /// Return true if all entries of the array are the same.
-   inline bool IsConstant() const;
+   bool IsConstant() const;
 
    /// Fill the entries of the array with the cumulative sum of the entries.
-   inline void PartialSum();
+   void PartialSum();
 
    /// Replace each entry of the array with its absolute value.
-   inline void Abs();
+   void Abs();
 
    /// Return the sum of all the array entries using the '+'' operator for class 'T'.
-   inline T Sum() const;
+   T Sum() const;
 
    /// Set all entries of the array to the provided constant.
    inline void operator=(const T &a);
@@ -803,14 +797,8 @@ template <typename T> template <typename CT>
 inline Array<T> &Array<T>::operator=(const Array<CT> &src)
 {
    SetSize(src.Size());
-
-   const bool use_dev = UseDevice() || src.UseDevice();
-   const auto x = src.Read(use_dev);
-   auto y = Write(use_dev);
-   mfem::forall_switch(use_dev, size, [=] MFEM_HOST_DEVICE (int i)
-   {
-      y[i] = x[i];
-   });
+   for (int i = 0; i < size; i++) { (*this)[i] = T(src[i]); }
+   return *this;
 }
 
 template <class T>
@@ -1026,24 +1014,19 @@ template <class T>
 inline void Array<T>::GetSubArray(int offset, int sa_size, Array<T> &sa) const
 {
    sa.SetSize(sa_size);
-   const bool use_dev = UseDevice() || sa.UseDevice();
-   const auto x = Read(use_dev);
-   auto y = sa.Write(use_dev);
-   mfem::forall_switch(use_dev, sa_size, [=] MFEM_HOST_DEVICE (int i)
+   for (int i = 0; i < sa_size; i++)
    {
-      y[i] = x[offset + i];
-   });
+      sa[i] = (*this)[offset+i];
+   }
 }
 
 template <class T>
 inline void Array<T>::operator=(const T &a)
 {
-   const bool use_dev = UseDevice();
-   auto x = Write(use_dev);
-   mfem::forall_switch(use_dev, size, [=] MFEM_HOST_DEVICE (int i)
+   for (int i = 0; i < size; i++)
    {
-      x[i] = a;
-   });
+      data[i] = a;
+   }
 }
 
 template <class T>
@@ -1052,153 +1035,6 @@ inline void Array<T>::Assign(const T *p)
    data.CopyFromHost(p, Size());
 }
 
-template <class T>
-inline void Array<T>::Print(std::ostream &os, int width) const
-{
-   for (int i = 0; i < size; i++)
-   {
-      os << data[i];
-      if ( !((i+1) % width) || i+1 == size )
-      {
-         os << '\n';
-      }
-      else
-      {
-         os << " ";
-      }
-   }
-}
-
-template <class T>
-inline void Array<T>::Save(std::ostream &os, int fmt) const
-{
-   if (fmt == 0)
-   {
-      os << size << '\n';
-   }
-   for (int i = 0; i < size; i++)
-   {
-      os << operator[](i) << '\n';
-   }
-}
-
-template <class T>
-void Array<T>::Load(std::istream &in, int fmt)
-{
-   if (fmt == 0)
-   {
-      int new_size;
-      in >> new_size;
-      SetSize(new_size);
-   }
-   for (int i = 0; i < size; i++)
-   {
-      in >> operator[](i);
-   }
-}
-
-template <class T>
-inline T Array<T>::Max() const
-{
-   MFEM_ASSERT(size > 0, "Array is empty with size " << size);
-
-   T max = operator[](0);
-   for (int i = 1; i < size; i++)
-   {
-      if (max < operator[](i))
-      {
-         max = operator[](i);
-      }
-   }
-
-   return max;
-}
-
-template <class T>
-inline T Array<T>::Min() const
-{
-   MFEM_ASSERT(size > 0, "Array is empty with size " << size);
-
-   T min = operator[](0);
-   for (int i = 1; i < size; i++)
-   {
-      if (operator[](i) < min)
-      {
-         min = operator[](i);
-      }
-   }
-
-   return min;
-}
-
-// Partial Sum
-template <class T>
-inline void Array<T>::PartialSum()
-{
-   auto data_ptr = ReadWrite(UseDevice());
-   InclusiveScan(UseDevice(), data_ptr, data_ptr, size);
-}
-
-template <class T>
-inline void Array<T>::Abs()
-{
-   static_assert(std::is_arithmetic<T>::value, "Use with arithmetic types!");
-   const bool useDevice = UseDevice();
-   const int N = size;
-   auto y = ReadWrite(useDevice);
-   mfem::forall_switch(useDevice, N, [=] MFEM_HOST_DEVICE (int i)
-   {
-      y[i] = std::abs(y[i]);
-   });
-}
-
-// Sum
-template <class T>
-inline T Array<T>::Sum() const
-{
-   T sum = static_cast<T>(0);
-   if (size > 0)
-   {
-      const auto m_data = Read(UseDevice());
-      reduce(size, sum, [=] MFEM_HOST_DEVICE(int i, T &r) { r += m_data[i]; },
-      /*  */ SumReducer<T> {}, UseDevice());
-   }
-   return sum;
-}
-
-template <class T>
-inline int Array<T>::IsSorted() const
-{
-   T val_prev = operator[](0), val;
-   for (int i = 1; i < size; i++)
-   {
-      val=operator[](i);
-      if (val < val_prev)
-      {
-         return 0;
-      }
-      val_prev = val;
-   }
-
-   return 1;
-}
-
-template <class T>
-inline bool Array<T>::IsConstant() const
-{
-   if (size < 2) { return true; }
-   const T v0 = data[0];
-   for (int i = 1; i < size; i++)
-   {
-      if (data[i] != v0)
-      {
-         return false;
-      }
-   }
-
-   return true;
-}
-
 
 template <class T>
 inline const T &Array2D<T>::operator()(int i, int j) const
@@ -1238,40 +1074,6 @@ inline T *Array2D<T>::operator[](int i)
    return &array1d[i*N];
 }
 
-template <class T>
-void Array2D<T>::Load(const char *filename, int fmt)
-{
-   std::ifstream in;
-   in.open(filename, std::ifstream::in);
-   MFEM_VERIFY(in.is_open(), "File " << filename << " does not exist.");
-   Load(in, fmt);
-   in.close();
-}
-
-template <class T>
-void Array2D<T>::Print(std::ostream &os, int width_)
-{
-   int height = this->NumRows();
-   int width  = this->NumCols();
-
-   for (int i = 0; i < height; i++)
-   {
-      os << "[row " << i << "]\n";
-      for (int j = 0; j < width; j++)
-      {
-         os << (*this)(i,j);
-         if ( (j+1) == width_ || (j+1) % width_ == 0 )
-         {
-            os << '\n';
-         }
-         else
-         {
-            os << ' ';
-         }
-      }
-   }
-}
-
 
 template <class T>
 inline void Swap(Array2D<T> &a, Array2D<T> &b)

general/reducers.hpp

@@ -12,6 +12,7 @@
 #ifndef MFEM_REDUCERS_HPP
 #define MFEM_REDUCERS_HPP
 
+#include "array.hpp"
 #include "forall.hpp"
 
 #include <cmath>
@@ -513,33 +514,6 @@ template<class B, class R> struct reduction_kernel
       }
    }
 };
-
-template <class T>
-class ReductionWorkspace
-{
-   Memory<T> workspace;
-
-   static ReductionWorkspace &Instance()
-   {
-      static ReductionWorkspace instance;
-      return instance;
-   }
-
-   ~ReductionWorkspace() { workspace.Delete(); }
-
-public:
-   static T *Get(int num_blocks)
-   {
-      ReductionWorkspace &instance = Instance();
-      if (instance.workspace.Capacity() < num_blocks)
-      {
-         instance.workspace.Delete();
-         instance.workspace.New(num_blocks, MemoryType::HOST_PINNED);
-      }
-      return instance.workspace;
-   }
-};
-
 }
 
 /**
@@ -555,7 +529,8 @@ public:
  @tparam T value_type to operate on
  */
 template <class T, class B, class R>
-void reduce(int N, T &res, B &&body, const R &reducer, bool use_dev)
+void reduce(int N, T &res, B &&body, const R &reducer, bool use_dev,
+            Array<T> &workspace)
 {
    if (N == 0)
    {
@@ -592,7 +567,13 @@ void reduce(int N, T &res, B &&body, const R &reducer, bool use_dev)
 
       red_type red{nullptr, std::forward<B>(body), reducer, N, items_per_thread};
       // allocate res to fit block_size entries
-      auto work = internal::ReductionWorkspace<T>::Get(nblocks);
+      auto mt = workspace.GetMemory().GetMemoryType();
+      if (mt != MemoryType::HOST_PINNED && mt != MemoryType::MANAGED)
+      {
+         mt = MemoryType::HOST_PINNED;
+      }
+      workspace.SetSize(nblocks, mt);
+      auto work = workspace.HostWrite();
       red.work = work;
       forall_2D(nblocks, block_size, 1, std::move(red));
       // wait for results

general/scan.hpp

@@ -28,37 +28,8 @@
 
 namespace mfem
 {
-
-namespace internal
-{
-class ScanWorkspace
-{
-   Memory<std::byte> workspace;
-   static ScanWorkspace &Instance()
-   {
-      static ScanWorkspace instance;
-      return instance;
-   }
-   ~ScanWorkspace() { workspace.Delete(); }
-public:
-   static std::byte *Get(int num_bytes)
-   {
-      ScanWorkspace &instance = Instance();
-      if (Size() < num_bytes)
-      {
-         instance.workspace.Delete();
-         instance.workspace.New(num_bytes);
-      }
-      return instance.workspace.Write(MemoryClass::DEVICE, Size());
-   }
-   static int Size()
-   {
-      return Instance().workspace.Capacity();
-   }
-};
-}
-
-/// Equivalent to InclusiveScan(use_dev, d_in, d_out, num_items, std::plus<>{})
+/// Equivalent to InclusiveScan(use_dev, d_in, d_out, num_items, workspace,
+/// std::plus<>{})
 template <class InputIt, class OutputIt>
 void InclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items)
 {
@@ -66,12 +37,12 @@ void InclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items)
 #if defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP)
    if (use_dev && mfem::Device::Allows(Backend::CUDA_MASK | Backend::HIP_MASK))
    {
-      using internal::ScanWorkspace;
-      size_t bytes = ScanWorkspace::Size();
-      if (bytes > 0)
+      static Array<std::byte> workspace;
+      size_t bytes = workspace.Size();
+      if (bytes)
       {
          auto err = MFEM_CUB_NAMESPACE::DeviceScan::InclusiveSum(
-                       ScanWorkspace::Get(bytes), bytes, d_in, d_out, num_items);
+                       workspace.Write(), bytes, d_in, d_out, num_items);
 #if defined(MFEM_USE_CUDA)
          if (err == cudaSuccess)
          {
@@ -86,12 +57,11 @@ void InclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items)
       }
       // try allocating a larger buffer
       bytes = 0;
-      // get size of buffer
       MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::InclusiveSum(
                         nullptr, bytes, d_in, d_out, num_items));
-      // resize buffer (in ScanWorkspace::Get) and try again
+      workspace.SetSize(bytes);
       MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::InclusiveSum(
-                        ScanWorkspace::Get(bytes), bytes, d_in, d_out, num_items));
+                        workspace.Write(), bytes, d_in, d_out, num_items));
       return;
    }
 #endif
@@ -131,13 +101,12 @@ void InclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
 #if defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP)
    if (use_dev && mfem::Device::Allows(Backend::CUDA_MASK | Backend::HIP_MASK))
    {
-      using internal::ScanWorkspace;
-      size_t bytes = ScanWorkspace::Size();
-      if (bytes > 0)
+      static Array<std::byte> workspace;
+      size_t bytes = workspace.Size();
+      if (bytes)
       {
          auto err = MFEM_CUB_NAMESPACE::DeviceScan::InclusiveScan(
-                       ScanWorkspace::Get(bytes), bytes, d_in, d_out, scan_op,
-                       num_items);
+                       workspace.Write(), bytes, d_in, d_out, scan_op, num_items);
 #if defined(MFEM_USE_CUDA)
          if (err == cudaSuccess)
          {
@@ -154,9 +123,9 @@ void InclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
       bytes = 0;
       MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::InclusiveScan(
                         nullptr, bytes, d_in, d_out, scan_op, num_items));
+      workspace.SetSize(bytes);
       MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::InclusiveScan(
-                        ScanWorkspace::Get(bytes), bytes, d_in, d_out, scan_op,
-                        num_items));
+                        workspace.Write(), bytes, d_in, d_out, scan_op, num_items));
       return;
    }
 #endif
@@ -195,13 +164,13 @@ void ExclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
 #if defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP)
    if (use_dev && mfem::Device::Allows(Backend::CUDA_MASK | Backend::HIP_MASK))
    {
-      using internal::ScanWorkspace;
-      size_t bytes = ScanWorkspace::Size();
+      static Array<std::byte> workspace;
+      size_t bytes = workspace.Size();
       if (bytes)
       {
          auto err = MFEM_CUB_NAMESPACE::DeviceScan::ExclusiveScan(
-                       ScanWorkspace::Get(bytes), bytes, d_in, d_out, scan_op,
-                       init_value, num_items);
+                       workspace.Write(), bytes, d_in, d_out, scan_op, init_value,
+                       num_items);
 #if defined(MFEM_USE_CUDA)
          if (err == cudaSuccess)
          {
@@ -218,9 +187,10 @@ void ExclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
       bytes = 0;
       MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::ExclusiveScan(
                         nullptr, bytes, d_in, d_out, scan_op, init_value, num_items));
+      workspace.SetSize(bytes);
       MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::ExclusiveScan(
-                        ScanWorkspace::Get(bytes), bytes, d_in, d_out, scan_op,
-                        init_value, num_items));
+                        workspace.Write(), bytes, d_in, d_out, scan_op, init_value,
+                        num_items));
       return;
    }
 #endif
@@ -243,7 +213,7 @@ void ExclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
 }
 
 /// Equivalent to ExclusiveScan(use_dev, d_in, d_out, num_items, init_value,
-/// std::plus<>{})
+/// workspace, std::plus<>{})
 template <class InputIt, class OutputIt, class T>
 void ExclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
                    T init_value)

linalg/vector.cpp

@@ -92,6 +92,18 @@ struct LpReducer
    }
 };
 
+static Array<real_t>& vector_workspace()
+{
+   static Array<real_t> instance;
+   return instance;
+}
+
+static Array<DevicePair<real_t, real_t>> &Lpvector_workspace()
+{
+   static Array<DevicePair<real_t, real_t>> instance;
+   return instance;
+}
+
 Vector::Vector(const Vector &v)
 {
    const int s = v.Size();
@@ -979,7 +991,7 @@ real_t Vector::Norml2() const
          }
       }
    },
-   L2Reducer{}, UseDevice());
+   L2Reducer{}, UseDevice(), Lpvector_workspace());
    // final answer
    return res.second * sqrt(res.first);
 }
@@ -994,7 +1006,7 @@ real_t Vector::Normlinf() const
    {
       r = fmax(r, fabs(m_data[i]));
    },
-   MaxReducer<real_t> {}, UseDevice());
+   MaxReducer<real_t> {}, UseDevice(), vector_workspace());
    return res;
 }
 
@@ -1008,7 +1020,7 @@ real_t Vector::Norml1() const
    {
       r += fabs(m_data[i]);
    },
-   SumReducer<real_t> {}, UseDevice());
+   SumReducer<real_t> {}, UseDevice(), vector_workspace());
    return res;
 }
 
@@ -1051,7 +1063,7 @@ real_t Vector::Normlp(real_t p) const
             }
          }
       },
-      LpReducer{p}, UseDevice());
+      LpReducer{p}, UseDevice(), Lpvector_workspace());
       // final answer
       return res.second * pow(res.first, 1.0 / p);
    } // end if p < infinity()
@@ -1084,7 +1096,7 @@ real_t Vector::operator*(const Vector &v) const
       {
          r += m_data[i] * v_data[i];
       },
-      SumReducer<real_t> {}, use_dev);
+      SumReducer<real_t> {}, use_dev, vector_workspace());
       return res;
    };
 
@@ -1155,7 +1167,7 @@ real_t Vector::Min() const
       {
          r = fmin(r, m_data[i]);
       },
-      MinReducer<real_t> {}, use_dev);
+      MinReducer<real_t> {}, use_dev, vector_workspace());
       return res;
    };
 
@@ -1201,7 +1213,7 @@ real_t Vector::Max() const
       {
          r = fmax(r, m_data[i]);
       },
-      MaxReducer<real_t> {}, use_dev);
+      MaxReducer<real_t> {}, use_dev, vector_workspace());
       return res;
    };
 
@@ -1236,7 +1248,7 @@ real_t Vector::Sum() const
    {
       r += m_data[i];
    },
-   SumReducer<real_t> {}, UseDevice());
+   SumReducer<real_t> {}, UseDevice(), vector_workspace());
    return res;
 }
 

miniapps/tools/gf-conv.cpp

@@ -0,0 +1,577 @@
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+static int dim;
+
+enum SCA_TYPE {INVALID_SCA_TYPE = -1,
+               H1_TYPE = 0,
+               L2_TYPE,
+               L2I_TYPE,
+               NUM_SCA_TYPES
+              };
+enum CONV_TYPE {INVALID_CONV_TYPE = -1,
+                PROJECTION = 0,
+                INTERPOLATION_OP,
+                SOLVE,
+                SOLVE_W_DBC,
+                NUM_CONV_TYPES
+               };
+
+FiniteElementCollection * GetFECollection(SCA_TYPE type, int p);
+ParFiniteElementSpace * GetFESpace(SCA_TYPE type, ParMesh &pmesh,
+                                   FiniteElementCollection &fec);
+void parseFieldNames(const char * field_name_c_str,
+                     vector<string> &field_names);
+
+string GetTypeName(SCA_TYPE type);
+string GetTypeShortName(SCA_TYPE type);
+string GetConvTypeName(CONV_TYPE type);
+string GetConvTypeShortName(CONV_TYPE type);
+
+void Projection(const ParGridFunction &v0, ParGridFunction &v1);
+void InterpolationOp(const ParGridFunction &v0, ParGridFunction &v1);
+void LeastSquares(SCA_TYPE t0, const ParGridFunction &v0,
+                  SCA_TYPE t1, ParGridFunction &v1);
+void LeastSquaresBC(SCA_TYPE t0, const ParGridFunction &v0,
+                    SCA_TYPE t1, ParGridFunction &v1,
+                    Coefficient &c);
+
+int main(int argc, char *argv[])
+{
+#ifdef MFEM_USE_MPI
+   Mpi::Init();
+   if (!Mpi::Root()) { mfem::out.Disable(); mfem::err.Disable(); }
+   Hypre::Init();
+#endif
+
+   // Parse command-line options.
+   const char *coll_name = NULL;
+   int cycle = 0;
+
+   const char *field_name_c_str = "ALL";
+
+   Array<int> orders;
+   Array<int> types;
+   Array<int> conv_types;
+   bool static_cond = false;
+   bool pa = false;
+   const char *device_config = "cpu";
+   bool visualization = 1;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&coll_name, "-r", "--root-file",
+                  "Set the VisIt data collection root file prefix.", true);
+   args.AddOption(&cycle, "-c", "--cycle", "Set the cycle index to read.");
+   args.AddOption(&field_name_c_str, "-fn", "--field-names",
+                  "List of field names to get values from.");
+   args.AddOption(&orders, "-o", "--final-order",
+                  "Finite element orders for each final field "
+                  "(an array of integers for multiple fields).");
+   args.AddOption(&types, "-t", "--final-type",
+                  "Set the basis type for the final fields: "
+                  "0-H1, 1-L2, 2-L2I, -1 loop over all.");
+   args.AddOption(&conv_types, "-ct", "--conversion-type",
+                  "Set the conversion schemes: "
+                  "0-Projection, 1-Interpolation Op, 2-Least Squares, "
+                  "3-Least Squares with BC, -1 loop over all.");
+   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
+                  "--no-static-condensation", "Enable static condensation.");
+   args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
+                  "--no-partial-assembly", "Enable Partial Assembly.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(mfem::out);
+      return 1;
+   }
+   args.PrintOptions(mfem::out);
+
+#ifdef MFEM_USE_MPI
+   VisItDataCollection dc(MPI_COMM_WORLD, coll_name);
+#else
+   VisItDataCollection dc(coll_name);
+#endif
+
+   dc.Load(cycle);
+
+   if (dc.Error() != DataCollection::NO_ERROR)
+   {
+      mfem::out << "Error loading VisIt data collection: " << coll_name << endl;
+      return 1;
+   }
+
+   dim = dc.GetMesh()->Dimension();
+   int spaceDim = dc.GetMesh()->SpaceDimension();
+
+   mfem::out << endl;
+   mfem::out << "Collection Name:    " << dc.GetCollectionName() << endl;
+   mfem::out << "Manifold Dimension: " << dim << endl;
+   mfem::out << "Space Dimension:    " << spaceDim << endl;
+   mfem::out << "Cycle:              " << dc.GetCycle() << endl;
+   mfem::out << "Time:               " << dc.GetTime() << endl;
+   mfem::out << "Time Step:          " << dc.GetTimeStep() << endl;
+   mfem::out << endl;
+
+   typedef DataCollection::FieldMapType fields_t;
+   const fields_t &fields = dc.GetFieldMap();
+   // Print the names of all fields.
+   mfem::out << "fields: [ ";
+   for (fields_t::const_iterator it = fields.begin(); it != fields.end(); ++it)
+   {
+      if (it != fields.begin()) { mfem::out << ", "; }
+      mfem::out << it->first;
+   }
+   mfem::out << " ]" << endl;
+
+   // Parsing desired field names
+   vector<string> field_names;
+   parseFieldNames(field_name_c_str, field_names);
+
+   if (field_names.size() == 1)
+   {
+      if (field_names[0] == "ALL")
+      {
+         fields_t::const_iterator it = fields.begin();
+         field_names[0] = it->first; it++;
+         for ( ; it != fields.end(); ++it)
+         {
+            field_names.push_back(it->first);
+         }
+      }
+   }
+
+   if (orders.Size() < field_names.size())
+   {
+      int size = orders.Size();
+      int order = (size > 0) ? orders[0] : 1;
+
+      orders.SetSize(field_names.size());
+      for (int i=size; i < field_names.size(); i++)
+      {
+         orders[i] = order;
+      }
+   }
+
+   if (types.Size() < field_names.size())
+   {
+      int size = types.Size();
+      int type = (size > 0) ? types[0] : 0;
+
+      types.SetSize(field_names.size());
+      for (int i=size; i < field_names.size(); i++)
+      {
+         types[i] = type;
+      }
+   }
+
+   if (conv_types.Size() < field_names.size())
+   {
+      int size = conv_types.Size();
+      int type = (size > 0) ? conv_types[0] : 0;
+
+      conv_types.SetSize(field_names.size());
+      for (int i=size; i < field_names.size(); i++)
+      {
+         conv_types[i] = type;
+      }
+   }
+
+   // Print field names to be extracted
+   mfem::out << "Extracting fields: ";
+   for (int i=0; i < field_names.size(); i++)
+   {
+      mfem::out << " \"" << field_names[i] << "\"";
+   }
+   mfem::out << endl;
+
+#ifdef MFEM_USE_MPI
+   ParMesh *mesh = dynamic_cast<ParMesh*>(dc.GetMesh());
+#else
+   Mesh *mesh = dc.GetMesh();
+#endif
+   if (mesh == NULL)
+   {
+      mfem::out << "Problem with mesh\n";
+      return 1;
+   }
+
+   int Ww = 300, Wh = 220, Fw = 3, Fh = 23, Ws = 15;
+
+   // Loop over all requested fields.
+   for (int i=0; i < field_names.size(); i++)
+   {
+#ifdef MFEM_USE_MPI
+      ParGridFunction *x0 = dc.GetParField(field_names[i]);
+#else
+      GridFunction *x0 = dc.GetField(field_names[i]);
+#endif
+      if (x0 == NULL)
+      {
+         mfem::out << "Problem with x0 for field \"" << field_names[i] << "\"\n";
+         continue;
+      }
+
+      int t0 = 0;
+
+      // nn. Send the solution by socket to a GLVis server.
+      if (visualization)
+      {
+         ostringstream oss;
+         oss << field_names[i];
+         char vishost[] = "localhost";
+         int  visport   = 19916;
+         socketstream sol_sock0(vishost, visport);
+#ifdef MFEM_USE_MPI
+         sol_sock0 << "parallel " << mesh->GetNRanks() << ' '
+                   << mesh->GetMyRank() << '\n';
+#endif
+         sol_sock0.precision(8);
+         sol_sock0 << "solution\n" << *mesh << *x0
+                   << "window_title '" << oss.str() << "'"
+                   << "window_geometry "
+                   << Ws * (t0) << " " << Ws * (t0) << " "
+                   << (int)(1.5 * Ww) << " " << (int)(1.5 * Wh)
+                   << flush;
+      }
+
+      int t1 = types[i];
+      FiniteElementCollection *fec1 = GetFECollection((SCA_TYPE)t1, orders[i]);
+      ParFiniteElementSpace   *fes1 = GetFESpace((SCA_TYPE)t1, *mesh, *fec1);
+
+      ParGridFunction *y1 = new ParGridFunction(fes1);
+
+      mfem::out << GetTypeName((SCA_TYPE)t1) << "(" << orders[i] << ")"
+                << ":" << endl;
+
+      int c01 = conv_types[i];
+      string cmnt = "";
+
+      switch ((CONV_TYPE)c01)
+      {
+         case PROJECTION:
+            Projection(*x0, *y1);
+            break;
+         case INTERPOLATION_OP:
+            cmnt = (t0 == (int)H1_TYPE) || (t0 == t1) ?
+                   "(should match projection)" : "(not expected to succeed)";
+            InterpolationOp(*x0, *y1);
+            break;
+         case SOLVE:
+            LeastSquares((SCA_TYPE)t0, *x0, (SCA_TYPE)t1, *y1);
+            break;
+         default:
+            *y1 = 0.0;
+      }
+
+      {
+         ostringstream oss;
+         oss << field_names[i] << "_" << GetConvTypeShortName((CONV_TYPE)c01)
+             << "_" << GetTypeShortName((SCA_TYPE)t1) << "_o" << orders[i];
+
+         dc.RegisterField(oss.str(), y1);
+      }
+
+      if (visualization)
+      {
+         ostringstream oss;
+         oss << GetConvTypeShortName((CONV_TYPE)c01) << "--> "
+             << GetTypeName((SCA_TYPE)t1)<< "(" << orders[i] << ")";
+         char vishost[] = "localhost";
+         int  visport   = 19916;
+         socketstream sol_sock1(vishost, visport);
+#ifdef MFEM_USE_MPI
+         sol_sock1 << "parallel " << mesh->GetNRanks() << ' '
+                   << mesh->GetMyRank() << '\n';
+#endif
+         sol_sock1.precision(8);
+         sol_sock1 << "solution\n" << *mesh << y1
+                   << "window_title '" << oss.str() << "'"
+                   << "window_geometry "
+                   << (int)((Ww + Fw) * (1.5 + c01) +
+                            Ws * (t0))
+                   << " " << (Wh + Fh) * (t1) + Ws * (t0)
+                   << " " << Ww << " " << Wh
+                   << flush;
+      }
+      mfem::out << endl;
+
+      // delete fes1;
+      // delete fec1;
+   }
+   dc.Save();
+
+   return 0;
+}
+
+FiniteElementCollection * GetFECollection(SCA_TYPE type, int p)
+{
+   switch (type)
+   {
+      case H1_TYPE:
+         return new H1_FECollection(p, dim);
+      case L2_TYPE:
+         return new L2_FECollection(p-1, dim);
+      case L2I_TYPE:
+         return new L2_FECollection(p-1, dim, BasisType::GaussLegendre,
+                                    FiniteElement::INTEGRAL);
+      default:
+         return NULL;
+   }
+}
+
+ParFiniteElementSpace * GetFESpace(SCA_TYPE type,
+                                   ParMesh &pmesh,
+                                   FiniteElementCollection &fec)
+{
+   return new ParFiniteElementSpace(&pmesh, &fec);
+}
+
+string GetTypeName(SCA_TYPE type)
+{
+   switch (type)
+   {
+      case H1_TYPE:
+         return "    H1";
+      case L2_TYPE:
+         return "    L2";
+      case L2I_TYPE:
+         return "    L2I";
+      default:
+         return "--";
+   }
+}
+
+string GetTypeShortName(SCA_TYPE type)
+{
+   switch (type)
+   {
+      case H1_TYPE:
+         return "H1";
+      case L2_TYPE:
+         return "L2";
+      case L2I_TYPE:
+         return "L2I";
+      default:
+         return "--";
+   }
+}
+
+string GetConvTypeName(CONV_TYPE type)
+{
+   switch (type)
+   {
+      case PROJECTION:
+         return "Projection            ";
+      case INTERPOLATION_OP:
+         return "Interpolation Operator";
+      case SOLVE:
+         return "Least Squares         ";
+      case SOLVE_W_DBC:
+         return "Least Squares with BC ";
+      default:
+         return "--";
+   }
+}
+
+string GetConvTypeShortName(CONV_TYPE type)
+{
+   switch (type)
+   {
+      case PROJECTION:
+         return "Proj";
+      case INTERPOLATION_OP:
+         return "Interp";
+      case SOLVE:
+         return "LS";
+      case SOLVE_W_DBC:
+         return "LSwBC";
+      default:
+         return "--";
+   }
+}
+
+void parseFieldNames(const char * field_name_c_str, vector<string> &field_names)
+{
+   string field_name_str(field_name_c_str);
+   string field_name;
+
+   for (string::iterator it=field_name_str.begin();
+        it!=field_name_str.end(); it++)
+   {
+      if (*it == '\\')
+      {
+         it++;
+         field_name.push_back(*it);
+      }
+      else if (*it == ' ')
+      {
+         if (!field_name.empty())
+         {
+            field_names.push_back(field_name);
+         }
+         field_name.clear();
+      }
+      else if (it == field_name_str.end() - 1)
+      {
+         field_name.push_back(*it);
+         field_names.push_back(field_name);
+      }
+      else
+      {
+         field_name.push_back(*it);
+      }
+   }
+   if (field_names.size() == 0)
+   {
+      field_names.push_back("ALL");
+   }
+}
+
+/** Perform a naive projection from one scalar field to another.
+
+   This scheme simply evaluates v0 at the interpolation points of v1.
+
+   If v0 has reduced continuity compared to v1 this can produce
+   results that depend on the order in which the elements are
+   traversed.
+
+   Suitable conversions:
+   H1 -> L2
+   H1 -> DG (same as L2)
+ */
+void Projection(const ParGridFunction &v0, ParGridFunction &v1)
+{
+   GridFunctionCoefficient v0Coef(&v0);
+   v1.ProjectCoefficient(v0Coef);
+}
+
+/** In theory this interpolation scheme should be equivalent to projection.
+
+    Building an interpolastion matrix could lead to computational
+    efficiency compared to simple projection if the operator will be
+    used several times.
+
+    Unfortunately this is broken for several combinations of source
+    and target fields.
+*/
+void InterpolationOp(const ParGridFunction &v0, ParGridFunction &v1)
+{
+   ParDiscreteLinearOperator op(v0.ParFESpace(), v1.ParFESpace());
+   op.AddDomainInterpolator(new IdentityInterpolator);
+   op.Assemble();
+   op.Finalize();
+
+   op.Mult(v0, v1);
+}
+
+/** Compute a least-squares best fit using the target basis functions.
+
+    This scheme is more difficult to setup and more computationally
+    expensive but the results can be significantly better than simple
+    projections.
+*/
+void LeastSquares(SCA_TYPE t0, const ParGridFunction &v0,
+                  SCA_TYPE t1, ParGridFunction &v1)
+{
+   ParFiniteElementSpace *fes0, *fes1;
+   fes0 = v0.ParFESpace();
+   fes1 = v1.ParFESpace();
+
+   ParMixedBilinearForm op(fes0, fes1);
+   op.AddDomainIntegrator(new MassIntegrator);
+   op.Assemble();
+   op.Finalize();
+
+   ParLinearForm b(v1.ParFESpace());
+   op.Mult(v0, b);
+
+   ParBilinearForm m(v1.ParFESpace());
+   m.AddDomainIntegrator(new MassIntegrator);
+   m.Assemble();
+   m.Finalize();
+
+   HypreParMatrix * M = m.ParallelAssemble();
+
+   HypreDiagScale diag(*M);
+   HyprePCG pcg(*M);
+   pcg.SetPreconditioner(diag);
+   pcg.SetTol(1e-12);
+   pcg.SetMaxIter(1000);
+
+   Vector B, X;
+   b.ParallelAssemble(B);
+
+   X.SetSize(v1.ParFESpace()->TrueVSize()); X = 0.0;
+   pcg.Mult(B, X);
+   v1.Distribute(X);
+
+   delete M;
+}
+
+/** Compute a least-squares best fit with boundary conditions.
+
+    This scheme is virtually identical to the previous one but it
+    makes use of boundary values, when available, to improve the
+    accuracy.  This scheme can produce significantly better results
+    when the normal derivative of the field is large near the
+    boundary.  This is particularly true when the field is
+    under-resolved near the boundary.
+*/
+void LeastSquaresBC(SCA_TYPE t0, const ParGridFunction &v0,
+                    SCA_TYPE t1, ParGridFunction &v1,
+                    Coefficient &c)
+{
+   ParFiniteElementSpace *fes0, *fes1;
+   fes0 = v0.ParFESpace();
+   fes1 = v1.ParFESpace();
+
+   ParMixedBilinearForm op(fes0, fes1);
+   op.AddDomainIntegrator(new MassIntegrator);
+   op.Assemble();
+   op.Finalize();
+
+   ParLinearForm b(v1.ParFESpace());
+   op.Mult(v0, b);
+
+   ParBilinearForm m(v1.ParFESpace());
+   m.AddDomainIntegrator(new MassIntegrator);
+   m.Assemble();
+   m.Finalize();
+
+   Array<int> ess_bdr;
+   Array<int> ess_tdof_list;
+   if (v1.ParFESpace()->GetParMesh()->bdr_attributes.Size())
+   {
+      ess_bdr.SetSize(v1.ParFESpace()->GetParMesh()->bdr_attributes.Max());
+      ess_bdr = 1;
+      v1.ParFESpace()->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+
+   if (t1 == H1_TYPE)
+   {
+      v1.ProjectBdrCoefficient(c, ess_bdr);
+   }
+
+   OperatorPtr M;
+   Vector B, X;
+   m.FormLinearSystem(ess_tdof_list, v1, b, M, X, B);
+
+   HypreDiagScale diag(*M.As<HypreParMatrix>());
+   HyprePCG pcg(*M.As<HypreParMatrix>());
+   pcg.SetPreconditioner(diag);
+   pcg.SetTol(1e-12);
+   pcg.SetMaxIter(1000);
+
+   pcg.Mult(B, X);
+   v1.Distribute(X);
+}

tests/unit/general/test_reduction.cpp

@@ -22,6 +22,7 @@ using namespace mfem;
 
 TEST_CASE("Reduce Sum", "[Reduction],[GPU]")
 {
+   Array<int> workspace;
    Array<int> a(1000);
    a.HostReadWrite();
    for (int i = 0; i < a.Size(); ++i)
@@ -35,7 +36,7 @@ TEST_CASE("Reduce Sum", "[Reduction],[GPU]")
       int res = 0;
       mfem::reduce(
       a.Size(), res, [=] MFEM_HOST_DEVICE(int i, int &r) { r += dptr[i]; },
-      SumReducer<int> {}, use_dev);
+      SumReducer<int> {}, use_dev, workspace);
       // correct for even-length summations
       int expected = (AsConst(a)[0] + AsConst(a)[a.Size() - 1]) * a.Size() / 2;
       CAPTURE(use_dev);
@@ -45,6 +46,7 @@ TEST_CASE("Reduce Sum", "[Reduction],[GPU]")
 
 TEST_CASE("Reduce Mult", "[Reduction],[GPU]")
 {
+   Array<long long> workspace;
    Array<long long> a(64);
    a.HostReadWrite();
    for (int i = 0; i < a.Size(); ++i)
@@ -62,7 +64,7 @@ TEST_CASE("Reduce Mult", "[Reduction],[GPU]")
          mfem::reduce(
             a.Size(), res,
          [=] MFEM_HOST_DEVICE(int i, long long &r) { r *= dptr[i]; },
-         MultReducer<long long> {}, use_dev);
+         MultReducer<long long> {}, use_dev, workspace);
          long long expected = 0;
          CAPTURE(use_dev);
          REQUIRE(res == expected);
@@ -74,7 +76,7 @@ TEST_CASE("Reduce Mult", "[Reduction],[GPU]")
          mfem::reduce(
             a.Size(), res,
          [=] MFEM_HOST_DEVICE(int i, long long &r) { r *= dptr[i]; },
-         MultReducer<long long> {}, use_dev);
+         MultReducer<long long> {}, use_dev, workspace);
          long long expected = 21936950640377856;
          CAPTURE(use_dev);
          REQUIRE(res == expected);
@@ -84,6 +86,7 @@ TEST_CASE("Reduce Mult", "[Reduction],[GPU]")
 
 TEST_CASE("Reduce BAnd", "[Reduction],[GPU]")
 {
+   Array<unsigned> workspace;
    Array<unsigned> a(10);
    SECTION("{ Bit unset }")
    {
@@ -105,7 +108,7 @@ TEST_CASE("Reduce BAnd", "[Reduction],[GPU]")
          mfem::reduce(
             a.Size(), res,
          [=] MFEM_HOST_DEVICE(int i, unsigned &r) { r &= dptr[i]; },
-         BAndReducer<unsigned> {}, use_dev);
+         BAndReducer<unsigned> {}, use_dev, workspace);
          CAPTURE(use_dev);
          REQUIRE(res == ((~1u) & ~(1u << unset_bit)));
          REQUIRE((res & (1u << unset_bit)) == 0);
@@ -129,7 +132,7 @@ TEST_CASE("Reduce BAnd", "[Reduction],[GPU]")
          mfem::reduce(
             a.Size(), res,
          [=] MFEM_HOST_DEVICE(int i, unsigned &r) { r &= dptr[i]; },
-         BAndReducer<unsigned> {}, use_dev);
+         BAndReducer<unsigned> {}, use_dev, workspace);
          CAPTURE(use_dev);
          REQUIRE(res == (1u << set_bit));
       }
@@ -138,6 +141,7 @@ TEST_CASE("Reduce BAnd", "[Reduction],[GPU]")
 
 TEST_CASE("Reduce BOr", "[Reduction],[GPU]")
 {
+   Array<unsigned> workspace;
    Array<unsigned> a(0x210);
    a.HostReadWrite();
    for (int i = 0; i < a.Size(); ++i)
@@ -153,7 +157,7 @@ TEST_CASE("Reduce BOr", "[Reduction],[GPU]")
       mfem::reduce(
          a.Size(), res,
       [=] MFEM_HOST_DEVICE(int i, unsigned &r) { r |= dptr[i]; },
-      BOrReducer<unsigned> {}, use_dev);
+      BOrReducer<unsigned> {}, use_dev, workspace);
       CAPTURE(use_dev);
       REQUIRE(res == 0x3ffu);
    }
@@ -161,6 +165,7 @@ TEST_CASE("Reduce BOr", "[Reduction],[GPU]")
 
 TEST_CASE("Reduce Min", "[Reduction],[GPU]")
 {
+   Array<int> workspace;
    Array<int> a(1000);
    auto hptr = a.HostReadWrite();
    for (int i = 0; i < a.Size(); ++i)
@@ -185,7 +190,7 @@ TEST_CASE("Reduce Min", "[Reduction],[GPU]")
             r = dptr[i];
          }
       },
-      MinReducer<int> {}, use_dev);
+      MinReducer<int> {}, use_dev, workspace);
       CAPTURE(use_dev);
       REQUIRE(res == -10);
    }
@@ -193,6 +198,7 @@ TEST_CASE("Reduce Min", "[Reduction],[GPU]")
 
 TEST_CASE("Reduce Max", "[Reduction],[GPU]")
 {
+   Array<int> workspace;
    Array<int> a(1000);
    auto hptr = a.HostReadWrite();
    for (int i = 0; i < a.Size(); ++i)
@@ -217,7 +223,7 @@ TEST_CASE("Reduce Max", "[Reduction],[GPU]")
             r = dptr[i];
          }
       },
-      MaxReducer<int> {}, use_dev);
+      MaxReducer<int> {}, use_dev, workspace);
       CAPTURE(use_dev);
       REQUIRE(res == 999 - 10);
    }
@@ -225,6 +231,7 @@ TEST_CASE("Reduce Max", "[Reduction],[GPU]")
 
 TEST_CASE("Reduce MinMax", "[Reduction],[GPU]")
 {
+   Array<DevicePair<int, int>> workspace;
    Array<int> a(1000);
    auto hptr = a.HostReadWrite();
    for (int i = 0; i < a.Size(); ++i)
@@ -255,7 +262,7 @@ TEST_CASE("Reduce MinMax", "[Reduction],[GPU]")
             r.second = dptr[i];
          }
       },
-      MinMaxReducer<int> {}, use_dev);
+      MinMaxReducer<int> {}, use_dev, workspace);
       CAPTURE(use_dev);
       REQUIRE(res.first == -10);
       REQUIRE(res.second == a.Size() - 11);
@@ -264,6 +271,7 @@ TEST_CASE("Reduce MinMax", "[Reduction],[GPU]")
 
 TEST_CASE("Reduce ArgMin", "[Reduction],[GPU]")
 {
+   Array<DevicePair<double, int>> workspace;
    Array<double> a(1000);
    auto hptr = a.HostReadWrite();
    for (int i = 0; i < a.Size(); ++i)
@@ -289,7 +297,7 @@ TEST_CASE("Reduce ArgMin", "[Reduction],[GPU]")
             r.second = i;
          }
       },
-      ArgMinReducer<double, int> {}, use_dev);
+      ArgMinReducer<double, int> {}, use_dev, workspace);
       CAPTURE(use_dev);
       REQUIRE(res.first == -10);
       REQUIRE(res.second >= 0);
@@ -300,6 +308,7 @@ TEST_CASE("Reduce ArgMin", "[Reduction],[GPU]")
 
 TEST_CASE("Reduce ArgMax", "[Reduction],[GPU]")
 {
+   Array<DevicePair<double, int>> workspace;
    Array<double> a(1000);
 
    auto hptr = a.HostReadWrite();
@@ -328,7 +337,7 @@ TEST_CASE("Reduce ArgMax", "[Reduction],[GPU]")
             r.second = i;
          }
       },
-      ArgMaxReducer<double, int> {}, use_dev);
+      ArgMaxReducer<double, int> {}, use_dev, workspace);
       CAPTURE(use_dev);
       REQUIRE(res.first == a.Size() - 11);
       REQUIRE(res.second >= 0);
@@ -339,6 +348,7 @@ TEST_CASE("Reduce ArgMax", "[Reduction],[GPU]")
 
 TEST_CASE("Reduce ArgMinMax", "[Reduction],[GPU]")
 {
+   Array<MinMaxLocScalar<double, int>> workspace;
    Array<double> a(1000);
    auto hptr = a.HostReadWrite();
    for (int i = 0; i < a.Size(); ++i)
@@ -373,7 +383,7 @@ TEST_CASE("Reduce ArgMinMax", "[Reduction],[GPU]")
             r.max_loc = i;
          }
       },
-      ArgMinMaxReducer<double, int> {}, use_dev);
+      ArgMinMaxReducer<double, int> {}, use_dev, workspace);
       CAPTURE(use_dev);
       REQUIRE(res.min_val == -10);
       REQUIRE(res.min_loc >= 0);