Reducing the two nested GMRES solvers to just one. Also implemented a complex version of AMS.

(this commit also improves some documentation in auxiliary.hpp)

CHANGELOG

@@ -11,6 +11,9 @@
 Version 4.2.1 (development)
 ===========================
 
+- Added high-order matrix-free auxiliary Maxwell solver for H(curl) problems,
+  as described in Barker and Kolev 2020 (https://doi.org/10.1002/nla.2348).
+
 - Added matrix-free GPU-enabled implementations of GradientInterpolator and
   IdentityInterpolator.
 
@@ -35,12 +38,6 @@ Version 4.2.1 (development)
 - Implemented a filter method for the Navier miniapp to stabilize highly
   turbulent flows in direct numerical simulation.
 
-- Added support for reading high-order Lagrange meshes in VTK format. Arbitrary-
-  orders and all element types are supported. See the VTK blog for more info:
-  https://blog.kitware.com/wp-content/uploads/2018/09/Source_Issue_43.pdf
-
-- Added support for reading VTK meshes in XML format.
-
 - Added partial assembly and device support to Example 25/25p, with diagonal
   preconditioning.
 
@@ -58,11 +55,6 @@ Version 4.2.1 (development)
   deprecated EvalSymmetric in MatrixCoefficient. Added DiagonalMatrixCoefficient
   for clarity, which is a typedef of VectorCoefficient.
 
-- Added support for AMG preconditioners for non-symmetric systems (e.g.
-  advection-dominated problems) using hypre's approximate ideal restriction
-  (AIR) AMG. Requires hypre version 2.14.0 or newer. Usage is illustrated in
-  example 9/9p.
-
 - Implemented an adaptive linear solver tolerance option for NewtonSolver based
   on the algorithm of Eisenstat and Walker.
 

NOTICE

@@ -29,5 +29,4 @@ license files. These software products and their licenses are as follows:
 * Catch++ (tests/unit/catch.hpp) -- Boost 1.0 license
 * Gecko (general/gecko.{cpp,hpp}) -- BSD 3-clause license
 * Picojson (fem/picojson.h) -- Custom 2-clause license
-* TinyXML2 (general/tinyxml2.{cpp,h}) -- zlib license
 * Zstr (general/zstr.hpp) -- MIT license

doc/CodeDocumentation.conf.in

@@ -819,9 +819,7 @@ RECURSIVE              = NO
 # run.
 
 EXCLUDE                = @MFEM_SOURCE_DIR@/config/_config.hpp \
-                         @MFEM_SOURCE_DIR@/config/get_hypre_version.cpp \
-                         @MFEM_SOURCE_DIR@/general/tinyxml2.h \
-                         @MFEM_SOURCE_DIR@/general/tinyxml2.cpp
+                         @MFEM_SOURCE_DIR@/config/get_hypre_version.cpp
 
 # The EXCLUDE_SYMLINKS tag can be used to select whether or not files or
 # directories that are symbolic links (a Unix file system feature) are excluded

examples/CMakeLists.txt

@@ -34,7 +34,6 @@ list(APPEND ALL_EXE_SRCS
   ex25.cpp
   ex26.cpp
   ex27.cpp
-  ex91.cpp
   )
 
 if (MFEM_USE_MPI)

examples/ex1p_complex.cpp

@@ -0,0 +1,435 @@
+//                       MFEM Example 3 - Parallel Version
+//
+// Compile with: make ex3p_complex
+//
+
+#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.
+double E_exact(const Vector &);
+void gradE_exact(const Vector &, Vector &);
+double f_exact(const Vector &);
+double freq = 1.0, kappa;
+int dim;
+
+
+#define COMPLEX_VERSION
+#define NEUMANN
+
+const double omega = 1.4;
+const double eps = 1.0e-8;
+
+
+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 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(&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(&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 (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      MPI_Finalize();
+      return 1;
+   }
+   if (myid == 0)
+   {
+      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 (myid == 0) { 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();
+   int sdim = mesh->SpaceDimension();
+
+   // 5. 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(1000./mesh->GetNE())/log(2.)/dim);
+      for (int l = 0; l < ref_levels; l++)
+      {
+         mesh->UniformRefinement();
+      }
+   }
+
+   // 6. 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 = 0;
+      for (int l = 0; l < par_ref_levels; l++)
+      {
+         pmesh->UniformRefinement();
+      }
+   }
+   pmesh->ReorientTetMesh();
+
+   // 7. Define a parallel finite element space on the parallel mesh. Here we
+   //    use the Nedelec finite elements of the specified order.
+   FiniteElementCollection *fec = new H1_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;
+   }
+
+   // 8. 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;
+
+#ifndef NEUMANN
+   if (pmesh->bdr_attributes.Size())
+   {
+      Array<int> ess_bdr(pmesh->bdr_attributes.Max());
+      ess_bdr = 1;
+      fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+#endif
+
+   //const double imscale = 0.0;
+   const double imscale = -omega;
+
+   Coefficient *im = new ConstantCoefficient(imscale);  // im part
+   //Coefficient *im = new ConstantCoefficient(0.0);  // im part
+
+   FunctionCoefficient E_coef(E_exact);
+   VectorFunctionCoefficient grad_E(sdim, gradE_exact);
+   ProductCoefficient omegaE(imscale, E_coef);
+
+   // 9. 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.
+   FunctionCoefficient f(f_exact);
+#ifdef COMPLEX_VERSION
+   ParComplexLinearForm *b = new ParComplexLinearForm(fespace);
+   b->AddDomainIntegrator(new DomainLFIntegrator(f), NULL);
+#ifdef NEUMANN
+   b->AddBoundaryIntegrator(NULL, new BoundaryNormalLFIntegrator(grad_E));
+#endif
+   b->AddBoundaryIntegrator(NULL, new BoundaryLFIntegrator(omegaE)); // im part
+#endif
+
+   b->Assemble();
+
+   // 10. 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);
+   */
+
+#ifdef COMPLEX_VERSION
+   // Complex version
+   ParComplexGridFunction x(fespace);
+   x = 0.0;
+   ConstantCoefficient E_im(0.0);
+   //x.ProjectBdrCoefficientTangent(E_Re, E_Im, ess_bdr);
+   x.ProjectCoefficient(E_coef, E_im);
+#endif
+
+   // 11. 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 *epscoef = new ConstantCoefficient(eps);
+   Coefficient *imabs = new ConstantCoefficient(fabs(imscale));  // im part
+
+#ifdef COMPLEX_VERSION
+   // Complex version
+   ParSesquilinearForm *a = new ParSesquilinearForm(fespace);
+   if (pa) { a->SetAssemblyLevel(AssemblyLevel::PARTIAL); }
+   a->AddDomainIntegrator(new DiffusionIntegrator(*muinv), NULL);
+   a->AddDomainIntegrator(new MassIntegrator(*epscoef), NULL);
+   a->AddBoundaryIntegrator(NULL, new MassIntegrator(*im));  // im part
+#endif
+
+   // 12. 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();
+
+   OperatorPtr A;
+   Vector B, X;
+   a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);
+
+   ParBilinearForm a_Re(fespace);
+   a_Re.AddDomainIntegrator(new DiffusionIntegrator(*muinv));
+   a_Re.AddDomainIntegrator(new MassIntegrator(*epscoef));
+
+   if (pa) { a_Re.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
+   a_Re.Assemble();
+
+   OperatorPtr A_Re;
+   a_Re.FormSystemMatrix(ess_tdof_list, A_Re);
+
+   ParBilinearForm a_Im(fespace);
+   a_Im.AddBoundaryIntegrator(new MassIntegrator(*imabs));
+   a_Im.Assemble();
+
+   OperatorPtr A_Im;
+   a_Im.FormSystemMatrix(ess_tdof_list, A_Im);
+
+   // 13. Solve the system AX=B using PCG with the AMS preconditioner from hypre
+   //     (in the full assembly case) or CG with Jacobi preconditioner (in the
+   //     partial assembly case).
+
+   Array<int> offsets(3);
+   offsets[0] = 0;
+   offsets[1] = fespace->GetTrueVSize();
+   offsets[2] = fespace->GetTrueVSize();
+   offsets.PartialSum();
+
+   //OperatorJacobiSmoother massJacobi(a_Im, ess_tdof_list);
+
+   StopWatch sw;
+   sw.Clear();
+   sw.Start();
+
+   if (pa) // Jacobi preconditioning in partial assembly mode
+   {
+      MFEM_VERIFY(false, "TODO");
+      //OperatorJacobiSmoother Jacobi(*a, ess_tdof_list);
+
+      CGSolver cg(MPI_COMM_WORLD);
+      cg.SetRelTol(1e-12);
+      cg.SetMaxIter(1000);
+      cg.SetPrintLevel(1);
+      cg.SetOperator(*A);
+      //cg.SetPreconditioner(Jacobi);
+      cg.Mult(B, X);
+   }
+   else
+   {
+      if (myid == 0)
+      {
+         cout << "Size of linear system: "
+              << A.As<HypreParMatrix>()->GetGlobalNumRows() << endl;
+      }
+
+      HypreBoomerAMG amg(*A_Re.As<HypreParMatrix>());
+
+#ifdef COMPLEX_VERSION
+      BlockDiagonalPreconditioner BlockDP(offsets);
+      BlockDP.SetDiagonalBlock(0, &amg);
+      BlockDP.SetDiagonalBlock(1, &amg);
+
+      Complex_PMHSS PMHSS(A_Re.Ptr(), A_Im.Ptr(), &BlockDP, NULL, 1.0);
+
+      ComplexOperator AspdComplex(A_Re.Ptr(), A_Im.Ptr(), false, false);
+
+      GMRESSolver PMHSSgmres(MPI_COMM_WORLD);
+      PMHSSgmres.SetPrintLevel(1);
+      PMHSSgmres.SetKDim(100);
+      PMHSSgmres.SetMaxIter(100);
+      PMHSSgmres.SetRelTol(1e-6);
+      PMHSSgmres.SetAbsTol(0.0);
+      PMHSSgmres.SetOperator(AspdComplex);
+      PMHSSgmres.SetPreconditioner(PMHSS);
+
+      GMRESSolver gmres(MPI_COMM_WORLD);
+      gmres.SetPrintLevel(1);
+      gmres.SetKDim(1000);
+      gmres.SetMaxIter(100);
+      gmres.SetRelTol(1e-8);
+      gmres.SetAbsTol(0.0);
+      gmres.SetOperator(*A);
+      //gmres.SetPreconditioner(BlockDP);
+      gmres.SetPreconditioner(PMHSS);
+      //gmres.SetPreconditioner(PMHSSgmres);
+#else
+      GMRESSolver gmres(MPI_COMM_WORLD);
+      gmres.SetPrintLevel(1);
+      gmres.SetKDim(1000);
+      gmres.SetMaxIter(100);
+      gmres.SetRelTol(1e-8);
+      gmres.SetAbsTol(0.0);
+      gmres.SetOperator(*A);
+      gmres.SetPreconditioner(ams);
+#endif
+
+      gmres.Mult(B, X);
+   }
+
+   sw.Stop();
+   mfem::out << "Total solve time " <<sw.RealTime() << endl;
+
+   // 14. Recover the parallel grid function corresponding to X. This is the
+   //     local finite element solution on each processor.
+   a->RecoverFEMSolution(X, *b, x);
+
+   // 15. Compute and print the L^2 norm of the error.
+   {
+#ifdef COMPLEX_VERSION
+      double err = x.real().ComputeL2Error(E_coef);
+#else
+      double err = x.ComputeL2Error(E_coef);
+#endif
+      if (myid == 0)
+      {
+         cout << "\n|| E_h - E ||_{L^2} = " << err << '\n' << endl;
+      }
+   }
+
+   // 16. 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);
+#ifdef COMPLEX_VERSION
+      x.real().Save(sol_ofs);
+#else
+      x.Save(sol_ofs);
+#endif
+   }
+
+   // 17. 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);
+#ifdef COMPLEX_VERSION
+      sol_sock << "solution\n" << *pmesh << x.real() << flush;
+#else
+      sol_sock << "solution\n" << *pmesh << x << flush;
+#endif
+   }
+
+   // 18. Free the used memory.
+   delete a;
+   delete muinv;
+   delete b;
+   delete fespace;
+   delete fec;
+   delete pmesh;
+
+   MPI_Finalize();
+
+   return 0;
+}
+
+#define VERSION_COS
+
+double E_exact(const Vector &x)
+{
+   if (dim == 3)
+   {
+#ifdef VERSION_COS
+      return cos(kappa * x(0)) * cos(kappa * x(1)) * cos(kappa * x(2));
+#else
+      return sin(kappa * x(0)) * sin(kappa * x(1)) * sin(kappa * x(2));
+#endif
+   }
+   else
+   {
+      return 0.0;
+   }
+}
+
+void gradE_exact(const Vector &x, Vector &grad)
+{
+   if (dim == 3)
+   {
+#ifdef VERSION_COS
+      grad(0) = -kappa * sin(kappa * x(0)) * cos(kappa * x(1)) * cos(kappa * x(2));
+      grad(1) = -kappa * sin(kappa * x(1)) * cos(kappa * x(0)) * cos(kappa * x(2));
+      grad(2) = -kappa * sin(kappa * x(2)) * cos(kappa * x(0)) * cos(kappa * x(1));
+#else
+      grad(0) = kappa * cos(kappa * x(0)) * sin(kappa * x(1)) * sin(kappa * x(2));
+      grad(1) = kappa * cos(kappa * x(1)) * sin(kappa * x(0)) * sin(kappa * x(2));
+      grad(2) = kappa * cos(kappa * x(2)) * sin(kappa * x(0)) * sin(kappa * x(1));
+#endif
+   }
+   else
+   {
+      MFEM_VERIFY(false, "");
+   }
+}
+
+// (grad u, grad v) + eps (u, v) = <grad u . n, v> - (div grad u, v) + eps (u, v)
+double f_exact(const Vector &x)
+{
+   if (dim == 3)
+   {
+      const double c = 3.0 * kappa * kappa;
+#ifdef VERSION_COS
+      return (eps + c) * cos(kappa * x(0)) * cos(kappa * x(1)) * cos(kappa * x(2));
+#else
+      return (eps + c) * sin(kappa * x(0)) * sin(kappa * x(1)) * sin(kappa * x(2));
+#endif
+   }
+   else
+   {
+      return 0.0;
+   }
+}

examples/ex3p.cpp

@@ -69,7 +69,10 @@ int main(int argc, char *argv[])
    bool static_cond = false;
    bool pa = false;
    const char *device_config = "cpu";
-   bool visualization = 1;
+   bool visualization = true;
+#ifdef MFEM_USE_AMGX
+   bool useAmgX = false;
+#endif
 
    OptionsParser args(argc, argv);
    args.AddOption(&mesh_file, "-m", "--mesh",
@@ -87,6 +90,11 @@ int main(int argc, char *argv[])
    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                   "--no-visualization",
                   "Enable or disable GLVis visualization.");
+#ifdef MFEM_USE_AMGX
+   args.AddOption(&useAmgX, "-amgx", "--useAmgX", "-no-amgx",
+                  "--no-useAmgX",
+                  "Enable or disable AmgX in MatrixFreeAMS.");
+#endif
 
    args.Parse();
    if (!args.Good())
@@ -159,9 +167,10 @@ int main(int argc, char *argv[])
    //    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;
+   Array<int> ess_bdr;
    if (pmesh->bdr_attributes.Size())
    {
-      Array<int> ess_bdr(pmesh->bdr_attributes.Max());
+      ess_bdr.SetSize(pmesh->bdr_attributes.Max());
       ess_bdr = 1;
       fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
    }
@@ -205,20 +214,20 @@ int main(int argc, char *argv[])
    Vector B, X;
    a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);
 
-   // 13. Solve the system AX=B using PCG with the AMS preconditioner from hypre
-   //     (in the full assembly case) or CG with Jacobi preconditioner (in the
-   //     partial assembly case).
-
-   if (pa) // Jacobi preconditioning in partial assembly mode
+   // 13. Solve the system AX=B using PCG with an AMS preconditioner.
+   if (pa)
    {
-      OperatorJacobiSmoother Jacobi(*a, ess_tdof_list);
-
+#ifdef MFEM_USE_AMGX
+      MatrixFreeAMS ams(*a, *A, *fespace, muinv, sigma, NULL, ess_bdr, useAmgX);
+#else
+      MatrixFreeAMS ams(*a, *A, *fespace, muinv, sigma, NULL, ess_bdr);
+#endif
       CGSolver cg(MPI_COMM_WORLD);
       cg.SetRelTol(1e-12);
       cg.SetMaxIter(1000);
       cg.SetPrintLevel(1);
       cg.SetOperator(*A);
-      cg.SetPreconditioner(Jacobi);
+      cg.SetPreconditioner(ams);
       cg.Mult(B, X);
    }
    else

examples/ex3p_complex.cpp

@@ -0,0 +1,517 @@
+//                       MFEM Example 3 - Parallel Version
+//
+// Compile with: make ex3p_complex
+//
+
+#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 curlE_exact(const Vector &, Vector &);
+void f_exact(const Vector &, Vector &);
+double freq = 1.0, kappa;
+int dim;
+
+
+#define COMPLEX_VERSION
+#define NEUMANN
+#define INDEFINITE
+
+const double omega = 1.4;
+
+
+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 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(&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(&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 (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      MPI_Finalize();
+      return 1;
+   }
+   if (myid == 0)
+   {
+      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 (myid == 0) { 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();
+   int sdim = mesh->SpaceDimension();
+
+   // 5. 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(1000./mesh->GetNE())/log(2.)/dim);
+      for (int l = 0; l < ref_levels; l++)
+      {
+         mesh->UniformRefinement();
+      }
+   }
+
+   // 6. 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 = 0;
+      for (int l = 0; l < par_ref_levels; l++)
+      {
+         pmesh->UniformRefinement();
+      }
+   }
+   pmesh->ReorientTetMesh();
+
+   // 7. 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;
+   }
+
+   // 8. 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;
+   Array<int> ess_bdr;
+   ess_bdr.SetSize(pmesh->bdr_attributes.Max());
+   ess_bdr = 0;
+
+#ifndef NEUMANN
+   if (pmesh->bdr_attributes.Size())
+   {
+      ess_bdr = 1;
+      fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+#endif
+
+   //const double imscale = 0.0;
+   const double imscale = omega;
+
+   Coefficient *im = new ConstantCoefficient(imscale);  // im part
+   //Coefficient *im = new ConstantCoefficient(0.0);  // im part
+
+   VectorFunctionCoefficient E_Re(sdim, E_exact);
+   VectorFunctionCoefficient curlE_Re(sdim, curlE_exact);
+
+   ScalarVectorProductCoefficient omegaE(imscale, E_Re);  // im part
+   //ScalarVectorProductCoefficient omegaE(0.0, E_Re);  // im part
+
+   // 9. 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);
+#ifdef COMPLEX_VERSION
+   ParComplexLinearForm *b = new ParComplexLinearForm(fespace);
+   b->AddDomainIntegrator(new VectorFEDomainLFIntegrator(f), NULL);
+   b->AddBoundaryIntegrator(NULL,
+                            new VectorFEDomainLFIntegrator(omegaE));  // im part
+#else
+   // Real version
+   ParLinearForm *b = new ParLinearForm(fespace);
+   b->AddDomainIntegrator(new VectorFEDomainLFIntegrator(f));
+#endif
+
+#ifdef NEUMANN
+   b->AddBoundaryIntegrator(new VectorFEBoundaryTangentLFIntegrator(curlE_Re),
+                            NULL);
+#endif
+
+   b->Assemble();
+
+   // 10. 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);
+   */
+
+#ifdef COMPLEX_VERSION
+   // Complex version
+   ParComplexGridFunction x(fespace);
+   x = 0.0;
+   Vector zero(sdim);
+   zero = 0.0;
+   VectorConstantCoefficient E_Im(zero);
+   //x.ProjectBdrCoefficientTangent(E_Re, E_Im, ess_bdr);
+   x.ProjectCoefficient(E_Re, E_Im);
+#else
+   ParGridFunction x(fespace);
+   x = 0.0;
+   x.ProjectCoefficient(E_Re);
+#endif
+
+   // 11. 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);
+#ifdef INDEFINITE
+   Coefficient *sigma = new ConstantCoefficient(
+      -omega*omega); // indefinite -, definite +
+#else
+   Coefficient *sigma = new ConstantCoefficient(
+      omega*omega); // indefinite -, definite +
+#endif
+   Coefficient *abssigma = new ConstantCoefficient(omega*omega);
+   Coefficient *imabs = new ConstantCoefficient(imscale);  // im part
+   //Coefficient *imabs = new ConstantCoefficient(0.0);  // im part
+   //Coefficient *im = new ConstantCoefficient(0.0);
+
+#ifdef COMPLEX_VERSION
+   // Complex version
+   ParSesquilinearForm *a = new ParSesquilinearForm(fespace);
+   if (pa) { a->SetAssemblyLevel(AssemblyLevel::PARTIAL); }
+   a->AddDomainIntegrator(new CurlCurlIntegrator(*muinv), NULL);
+   //a->AddDomainIntegrator(new VectorFEMassIntegrator(*sigma), new VectorFEMassIntegrator(*im));
+   a->AddDomainIntegrator(new VectorFEMassIntegrator(*sigma), NULL);
+   a->AddBoundaryIntegrator(NULL, new VectorFEMassIntegrator(*im));  // im part
+#else
+   // Real version
+   ParBilinearForm *a = new ParBilinearForm(fespace);
+   if (pa) { a->SetAssemblyLevel(AssemblyLevel::PARTIAL); }
+   a->AddDomainIntegrator(new CurlCurlIntegrator(*muinv));
+   //a->AddDomainIntegrator(new VectorFEMassIntegrator(*sigma), new VectorFEMassIntegrator(*im));
+   a->AddDomainIntegrator(new VectorFEMassIntegrator(*sigma));
+#endif
+
+   // 12. 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();
+
+   OperatorPtr A;
+   Vector B, X;
+   a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);
+
+   ParBilinearForm a_Re(fespace);
+   a_Re.AddDomainIntegrator(new CurlCurlIntegrator(*muinv));
+   a_Re.AddDomainIntegrator(new VectorFEMassIntegrator(*abssigma));
+
+   //if (pa) { a_Re.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
+   a_Re.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+   a_Re.Assemble();
+
+   OperatorPtr A_Re;
+   a_Re.FormSystemMatrix(ess_tdof_list, A_Re);
+
+   ParBilinearForm a_Im(fespace);
+   a_Im.AddBoundaryIntegrator(new VectorFEMassIntegrator(*imabs));
+   a_Im.Assemble();
+
+   OperatorPtr A_Im;
+   a_Im.FormSystemMatrix(ess_tdof_list, A_Im);
+
+   // 13. Solve the system AX=B using PCG with the AMS preconditioner from hypre
+   //     (in the full assembly case) or CG with Jacobi preconditioner (in the
+   //     partial assembly case).
+
+   Array<int> offsets(3);
+   offsets[0] = 0;
+   offsets[1] = fespace->GetTrueVSize();
+   offsets[2] = fespace->GetTrueVSize();
+   offsets.PartialSum();
+
+   //OperatorJacobiSmoother massJacobi(a_Im, ess_tdof_list);
+
+   StopWatch sw;
+   sw.Clear();
+   sw.Start();
+
+   if (pa) // Jacobi preconditioning in partial assembly mode
+   {
+      MFEM_VERIFY(false, "TODO");
+      //OperatorJacobiSmoother Jacobi(*a, ess_tdof_list);
+
+      CGSolver cg(MPI_COMM_WORLD);
+      cg.SetRelTol(1e-12);
+      cg.SetMaxIter(1000);
+      cg.SetPrintLevel(1);
+      cg.SetOperator(*A);
+      //cg.SetPreconditioner(Jacobi);
+      cg.Mult(B, X);
+   }
+   else
+   {
+      if (myid == 0)
+      {
+         cout << "Size of linear system: "
+              << A.As<HypreParMatrix>()->GetGlobalNumRows() << endl;
+      }
+
+      //HypreAMS ams(*A_Re.As<HypreParMatrix>(), fespace);
+
+      // One option is to use the standard real-valued MatrixFreeAMS to precondition
+      // the real part of the complex system in the PMHSS preconditioner (BlockDiagonalPreconditioner).
+      // Another option is to use complex MatrixFreeAMS to precondition the
+      // complex system without PMHSS and without a BlockDiagonalPreconditioner.
+      //#define COMPLEX_AMS
+
+#ifdef MFEM_USE_AMGX
+      bool useAmgX = false;
+      cout << "Built with AMGX, using AMGX " << useAmgX << endl;
+      MatrixFreeAMS ams(a_Re, *A_Re, *fespace, muinv, abssigma, im, imabs, NULL,
+                        ess_bdr, useAmgX);
+      MatrixFreeAMS ams(a_Re, *A_Re, *fespace, muinv, abssigma, NULL, NULL, ess_bdr,
+                        useAmgX);
+#ifdef COMPLEX_AMS
+      MFEM_VERIFY(false, "TODO");
+#endif
+
+#else
+      cout << "Not built with AMGX" << endl;
+#ifdef COMPLEX_AMS
+      MatrixFreeAMS ams(a_Re, *A_Re, A.Ptr(), *fespace, muinv, abssigma, im, imabs,
+                        NULL, ess_bdr);
+#else
+      MatrixFreeAMS ams(a_Re, *A_Re, NULL, *fespace, muinv, abssigma, NULL, NULL,
+                        NULL, ess_bdr);
+#endif
+#endif
+
+#ifdef COMPLEX_VERSION
+
+#ifdef COMPLEX_AMS
+      //MFEM_VERIFY(false, "TODO");
+#else
+      BlockDiagonalPreconditioner BlockDP(offsets);
+      BlockDP.SetDiagonalBlock(0, &ams);
+      BlockDP.SetDiagonalBlock(1, &ams);
+
+      /*
+      BlockDiagonalPreconditioner BlockDP_Im(offsets);
+      BlockDP_Im.SetDiagonalBlock(0, &massJacobi); // TODO: this won't work if it has zeros on diagonal
+      BlockDP_Im.SetDiagonalBlock(1, &massJacobi);
+      */
+
+      //Complex_PMHSS PMHSS(A_Re, A_Im, &BlockDP, &BlockDP_Im);
+      //Complex_PMHSS PMHSS(A_Re, A_Im, &BlockDP, NULL, 2.0 * omega);
+      //Complex_PMHSS PMHSS(A_Re, A_Im, &BlockDP, NULL, omega);
+      Complex_PMHSS PMHSS(A_Re.Ptr(), A_Im.Ptr(), &BlockDP, NULL, 1.0);
+
+      ComplexOperator AspdComplex(A_Re.Ptr(), A_Im.Ptr(), false, false);
+
+      GMRESSolver PMHSSgmres(MPI_COMM_WORLD);
+      PMHSSgmres.SetPrintLevel(1);
+      PMHSSgmres.SetKDim(100);
+      PMHSSgmres.SetMaxIter(100);
+      PMHSSgmres.SetRelTol(1e-6);
+      PMHSSgmres.SetAbsTol(0.0);
+      PMHSSgmres.SetOperator(AspdComplex);
+      PMHSSgmres.SetPreconditioner(PMHSS);
+#endif
+
+      GMRESSolver gmres(MPI_COMM_WORLD);
+      gmres.SetPrintLevel(1);
+      gmres.SetKDim(1000);
+      gmres.SetMaxIter(100);
+      gmres.SetRelTol(1e-8);
+      gmres.SetAbsTol(0.0);
+      gmres.SetOperator(*A);
+      //gmres.SetPreconditioner(BlockDP);
+#ifdef COMPLEX_AMS
+      //MFEM_VERIFY(false, "TODO");
+      gmres.SetPreconditioner(ams);
+#else
+      gmres.SetPreconditioner(PMHSS);
+      //gmres.SetPreconditioner(PMHSSgmres);
+#endif
+
+#else
+      GMRESSolver gmres(MPI_COMM_WORLD);
+      gmres.SetPrintLevel(1);
+      gmres.SetKDim(1000);
+      gmres.SetMaxIter(100);
+      gmres.SetRelTol(1e-8);
+      gmres.SetAbsTol(0.0);
+      gmres.SetOperator(*A);
+      gmres.SetPreconditioner(ams);
+#endif
+
+      gmres.Mult(B, X);
+   }
+
+   sw.Stop();
+   mfem::out << "Total solve time " <<sw.RealTime() << endl;
+
+   // 14. Recover the parallel grid function corresponding to X. This is the
+   //     local finite element solution on each processor.
+   a->RecoverFEMSolution(X, *b, x);
+
+   // 15. Compute and print the L^2 norm of the error.
+   {
+#ifdef COMPLEX_VERSION
+      double err = x.real().ComputeL2Error(E_Re);
+#else
+      double err = x.ComputeL2Error(E_Re);
+#endif
+      if (myid == 0)
+      {
+         cout << "\n|| E_h - E ||_{L^2} = " << err << '\n' << endl;
+      }
+   }
+
+   // 16. 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);
+#ifdef COMPLEX_VERSION
+      x.real().Save(sol_ofs);
+#else
+      x.Save(sol_ofs);
+#endif
+   }
+
+   // 17. 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);
+#ifdef COMPLEX_VERSION
+      sol_sock << "solution\n" << *pmesh << x.real() << flush;
+#else
+      sol_sock << "solution\n" << *pmesh << x << flush;
+#endif
+   }
+
+   // 18. 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 curlE_exact(const Vector &x, Vector &curl)
+{
+   if (dim == 3)
+   {
+      curl(0) = kappa * cos(kappa * x(2));
+      curl(1) = kappa * cos(kappa * x(0));
+      curl(2) = kappa * cos(kappa * x(1));
+   }
+   else
+   {
+      MFEM_VERIFY(false, "");
+   }
+}
+
+void f_exact(const Vector &x, Vector &f)
+{
+   if (dim == 3)
+   {
+      // indefinite -m, definite +m
+      const double c = kappa * kappa;
+#ifdef INDEFINITE
+      const double m = -omega * omega;
+#else
+      const double m = omega * omega;
+#endif
+      f(0) = (c + m) * sin(kappa * x(1));
+      f(1) = (c + m) * sin(kappa * x(2));
+      f(2) = (c + m) * 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/ex9.cpp

@@ -290,16 +290,15 @@ int main(int argc, char *argv[])
       k.SetAssemblyLevel(AssemblyLevel::FULL);
    }
    m.AddDomainIntegrator(new MassIntegrator);
-   constexpr double alpha = -1.0;
-   k.AddDomainIntegrator(new ConvectionIntegrator(velocity, alpha));
+   k.AddDomainIntegrator(new ConvectionIntegrator(velocity, -1.0));
    k.AddInteriorFaceIntegrator(
-      new NonconservativeDGTraceIntegrator(velocity, alpha));
+      new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
    k.AddBdrFaceIntegrator(
-      new NonconservativeDGTraceIntegrator(velocity, alpha));
+      new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
 
    LinearForm b(&fes);
    b.AddBdrFaceIntegrator(
-      new BoundaryFlowIntegrator(inflow, velocity, alpha));
+      new BoundaryFlowIntegrator(inflow, velocity, -1.0, -0.5));
 
    m.Assemble();
    int skip_zeros = 0;

examples/ex91.cpp

@@ -1,672 +0,0 @@
-#include"mfem.hpp"
-
-#include<memory>
-#include<iostream>
-#include<fstream>
-
-namespace mfem {
-
-
-class LinDiffQFunc
-{
-public:
-    LinDiffQFunc(mfem::Coefficient& dd, mfem::Coefficient& ll, double gg_,
-                 double pp0_, double pp1_):diff(dd),load(ll),gg(gg_),pp0(pp0_),pp1(pp1_)
-    {
-
-    }
-
-    double QEnergy(ElementTransformation &T,
-                  const IntegrationPoint &ip,
-                  mfem::Vector& param, mfem::Vector& uu)
-    {
-        double dd=diff.Eval(T,ip);
-        double ll=load.Eval(T,ip);
-
-        double rho0=param[0];
-        double rho1=param[1];
-        double fd=dd*std::pow(rho0,pp0)*std::pow(rho1,pp1);
-
-        double rez = 0.5*(uu[0]*uu[0]+uu[1]*uu[1]+uu[2]*uu[2])*fd
-                        + 0.5*gg*uu[3]*uu[3] -uu[3]*ll;
-
-        return rez;
-    }
-
-    void QResidual(ElementTransformation &T,
-                   const IntegrationPoint &ip,
-                   mfem::Vector& param, mfem::Vector& uu, mfem::Vector& rr)
-    {
-        rr.SetSize(4);
-        double dd=diff.Eval(T,ip);
-        double ll=load.Eval(T,ip);
-
-        double rho0=param[0];
-        double rho1=param[1];
-
-        double fd=dd*std::pow(rho0,pp0)*std::pow(rho1,pp1);
-
-        rr[0]=uu[0]*fd;
-        rr[1]=uu[1]*fd;
-        rr[2]=uu[2]*fd;
-        rr[3]=gg*uu[3]-ll;
-    }
-
-    void AQResidual(ElementTransformation &T,
-                    const IntegrationPoint &ip,
-                    mfem::Vector& param,
-                    mfem::Vector& uu, mfem::Vector& aa, mfem::Vector& rr)
-    {
-         rr.SetSize(2);
-         double dd=diff.Eval(T,ip);
-         double ll=load.Eval(T,ip);
-
-         double rho0=param[0];
-         double rho1=param[1];
-
-         double fd0=dd*pp0*std::pow(rho0,pp0-1.0)*std::pow(rho1,pp1);
-         double fd1=dd*std::pow(rho0,pp0)*pp1*std::pow(rho1,pp1-1.0);
-
-         rr[0] = (aa[0]*uu[0]+aa[1]*uu[1]+aa[2]*uu[2])*fd0;
-         rr[1] = (aa[0]*uu[0]+aa[1]*uu[1]+aa[2]*uu[2])*fd1;
-
-    }
-
-    void QGradResidual(ElementTransformation &T,
-                       const IntegrationPoint &ip,
-                       mfem::Vector& param, mfem::Vector& uu, mfem::DenseMatrix& hh)
-    {
-        hh.SetSize(4);
-        double dd=diff.Eval(T,ip);
-        //double ll=load.Eval(T,ip);
-
-
-        double rho0=param[0];
-        double rho1=param[1];
-
-        double fd=dd*std::pow(rho0,pp0)*std::pow(rho1,pp1);
-        hh=0.0;
-
-        hh(0,0)=fd;
-        hh(1,1)=fd;
-        hh(2,2)=fd;
-        hh(3,3)=gg;
-    }
-
-
-private:
-    mfem::Coefficient& diff;
-    mfem::Coefficient& load;
-    double gg;
-    double pp0;
-    double pp1;
-};
-
-
-class PrmBlockLSFEMDiffusion: public PrmBlockNonlinearFormIntegrator
-{
-public:
-    PrmBlockLSFEMDiffusion(LinDiffQFunc& qfun_)
-    {
-        qfunc=&qfun_;
-    }
-
-    /// Compute the local energy
-    virtual double GetElementEnergy(const Array<const FiniteElement *>&el,
-                                    const Array<const FiniteElement *>&pel,
-                                    ElementTransformation &Tr,
-                                    const Array<const Vector *>&elfun,
-                                    const Array<const Vector *>&pelfun)
-    {
-        int dof_u0 = el[0]->GetDof();
-        int dof_r0 = pel[0]->GetDof();
-        int dof_r1 = pel[1]->GetDof();
-
-        int dim = el[0]->GetDim();
-        int spaceDim = Tr.GetSpaceDim();
-        if (dim != spaceDim)
-        {
-            mfem::mfem_error(" PrmBlockLSFEMDiffusion::GetElementEnergy"
-                             " is not defined on manifold meshes");
-        }
-
-        //shape functions
-        Vector shu0(dof_u0);
-        Vector shr0(dof_r0);
-        Vector shr1(dof_r1);
-        DenseMatrix dsu0(dof_u0,dim);
-        DenseMatrix B(dof_u0, 4);
-        B=0.0;
-
-        double w;
-
-        Vector param(2); param=0.0;
-        Vector uu(4); uu=0.0;
-
-        double energy =0.0;
-
-        const IntegrationRule *ir = nullptr;
-        if(ir==nullptr){
-            int order= 2 * el[0]->GetOrder() + Tr.OrderGrad(el[0])
-                                +pel[0]->GetOrder()+pel[1]->GetOrder();
-            ir=&IntRules.Get(Tr.GetGeometryType(),order);
-        }
-
-        for (int i = 0; i < ir->GetNPoints(); i++)
-        {
-            const IntegrationPoint &ip = ir->IntPoint(i);
-            Tr.SetIntPoint(&ip);
-            w=Tr.Weight();
-            w = ip.weight * w;
-
-            el[0]->CalcPhysDShape(Tr,dsu0);
-            el[0]->CalcPhysShape(Tr,shu0);
-            pel[0]->CalcPhysShape(Tr,shr0);
-            pel[1]->CalcPhysShape(Tr,shr1);
-
-            param[0]=shr0*(*pelfun[0]);
-            param[1]=shr1*(*pelfun[1]);
-
-            //set the matrix B
-            for(int jj=0;jj<dim;jj++)
-            {
-                B.SetCol(jj,dsu0.GetColumn(jj));
-            }
-            B.SetCol(3,shu0);
-            B.MultTranspose(*elfun[0],uu);
-            energy=energy+w * qfunc->QEnergy(Tr,ip,param,uu);
-        }
-        return energy;
-    }
-
-    /// Perform the local action of the BlockNonlinearFormIntegrator
-    virtual void AssembleElementVector(const Array<const FiniteElement *> &el,
-                                       const Array<const FiniteElement *>&pel,
-                                       ElementTransformation &Tr,
-                                       const Array<const Vector *> &elfun,
-                                       const Array<const Vector *>&pelfun,
-                                       const Array<Vector *> &elvec)
-    {
-        int dof_u0 = el[0]->GetDof();
-        int dof_r0 = pel[0]->GetDof();
-        int dof_r1 = pel[1]->GetDof();
-
-        int dim = el[0]->GetDim();
-
-        elvec[0]->SetSize(dof_u0);
-        *elvec[0]=0.0;
-        int spaceDim = Tr.GetSpaceDim();
-        if (dim != spaceDim)
-        {
-            mfem::mfem_error(" PrmBlockLSFEMDiffusion::AssembleElementVector"
-                             " is not defined on manifold meshes");
-        }
-
-        //shape functions
-        Vector shu0(dof_u0);
-        Vector shr0(dof_r0);
-        Vector shr1(dof_r1);
-        DenseMatrix dsu0(dof_u0,dim);
-        DenseMatrix B(dof_u0, 4);
-        B=0.0;
-
-        double w;
-
-        Vector param(2); param=0.0;
-        Vector uu(4); uu=0.0;
-        Vector rr;
-        Vector lvec; lvec.SetSize(dof_u0);
-
-        const IntegrationRule *ir = nullptr;
-        if(ir==nullptr){
-            int order= 2 * el[0]->GetOrder() + Tr.OrderGrad(el[0])
-                                +pel[0]->GetOrder()+pel[1]->GetOrder();
-            ir=&IntRules.Get(Tr.GetGeometryType(),order);
-        }
-
-        for (int i = 0; i < ir->GetNPoints(); i++)
-        {
-            const IntegrationPoint &ip = ir->IntPoint(i);
-            Tr.SetIntPoint(&ip);
-            w=Tr.Weight();
-            w = ip.weight * w;
-
-            el[0]->CalcPhysDShape(Tr,dsu0);
-            el[0]->CalcPhysShape(Tr,shu0);
-            pel[0]->CalcPhysShape(Tr,shr0);
-            pel[1]->CalcPhysShape(Tr,shr1);
-
-            param[0]=shr0*(*pelfun[0]);
-            param[1]=shr1*(*pelfun[1]);
-
-            //set the matrix B
-            for(int jj=0;jj<dim;jj++)
-            {
-                B.SetCol(jj,dsu0.GetColumn(jj));
-            }
-            B.SetCol(3,shu0);
-            B.MultTranspose(*elfun[0],uu);
-            qfunc->QResidual(Tr,ip,param, uu, rr);
-
-            B.Mult(rr,lvec);
-            elvec[0]->Add(w,lvec);
-        }
-
-    }
-
-    virtual void AssembleFaceVector(const Array<const FiniteElement *> &el1,
-                                    const Array<const FiniteElement *> &el2,
-                                    const Array<const FiniteElement *> &pel1,
-                                    const Array<const FiniteElement *> &pel2,
-                                    FaceElementTransformations &Tr,
-                                    const Array<const Vector *> &elfun,
-                                    const Array<const Vector *>&pelfun,
-                                    const Array<Vector *> &elvect)
-    {
-
-    }
-
-    /// Assemble the local gradient matrix
-    virtual void AssembleElementGrad(const Array<const FiniteElement*> &el,
-                                     const Array<const FiniteElement *>&pel,
-                                     ElementTransformation &Tr,
-                                     const Array<const Vector *> &elfun,
-                                     const Array<const Vector *>&pelfun,
-                                     const Array2D<DenseMatrix *> &elmats)
-    {
-        int dof_u0 = el[0]->GetDof();
-        int dof_r0 = pel[0]->GetDof();
-        int dof_r1 = pel[1]->GetDof();
-
-        int dim = el[0]->GetDim();
-
-        //elmats[0]->Size(dof_u0, dof_u0);
-        //*elmats[0]=0.0;
-
-        DenseMatrix* K=elmats(0,0);
-        K->SetSize(dof_u0,dof_u0);
-        (*K)=0.0;
-
-        int spaceDim = Tr.GetSpaceDim();
-        if (dim != spaceDim)
-        {
-            mfem::mfem_error(" PrmBlockLSFEMDiffusion::AssembleElementVector"
-                             " is not defined on manifold meshes");
-        }
-
-        //shape functions
-        Vector shu0(dof_u0);
-        Vector shr0(dof_r0);
-        Vector shr1(dof_r1);
-        DenseMatrix dsu0(dof_u0,dim);
-        DenseMatrix B(dof_u0, 4);
-        DenseMatrix A(dof_u0, 4);
-        B=0.0;
-
-        double w;
-
-        Vector param(2); param=0.0;
-        Vector uu(4); uu=0.0;
-        DenseMatrix hh;
-        Vector lvec; lvec.SetSize(dof_u0);
-
-        const IntegrationRule *ir = nullptr;
-        if(ir==nullptr){
-            int order= 2 * el[0]->GetOrder() + Tr.OrderGrad(el[0])
-                                +pel[0]->GetOrder()+pel[1]->GetOrder();
-            ir=&IntRules.Get(Tr.GetGeometryType(),order);
-        }
-
-        for (int i = 0; i < ir->GetNPoints(); i++)
-        {
-            const IntegrationPoint &ip = ir->IntPoint(i);
-            Tr.SetIntPoint(&ip);
-            w = Tr.Weight();
-            w = ip.weight * w;
-
-            el[0]->CalcPhysDShape(Tr,dsu0);
-            el[0]->CalcPhysShape(Tr,shu0);
-            pel[0]->CalcPhysShape(Tr,shr0);
-            pel[1]->CalcPhysShape(Tr,shr1);
-
-            param[0]=shr0*(*pelfun[0]);
-            param[1]=shr1*(*pelfun[1]);
-
-            //set the matrix B
-            for(int jj=0;jj<dim;jj++)
-            {
-                B.SetCol(jj,dsu0.GetColumn(jj));
-            }
-            B.SetCol(3,shu0);
-            B.MultTranspose(*elfun[0],uu);
-            qfunc->QGradResidual(Tr,ip,param,uu,hh);
-            Mult(B,hh,A);
-            AddMult_a_ABt(w,A,B,*K);
-        }
-    }
-
-    virtual void AssembleFaceGrad(const Array<const FiniteElement *>&el1,
-                                  const Array<const FiniteElement *>&el2,
-                                  const Array<const FiniteElement *> &pel1,
-                                  const Array<const FiniteElement *> &pel2,
-                                  FaceElementTransformations &Tr,
-                                  const Array<const Vector *> &elfun,
-                                  const Array<const Vector *>&pelfun,
-                                  const Array2D<DenseMatrix *> &elmats)
-    {
-
-    }
-
-    virtual void AssemblePrmElementVector(const Array<const FiniteElement *> &el,
-                                          const Array<const FiniteElement *> &pel,
-                                          ElementTransformation &Tr,
-                                          const Array<const Vector *> &elfun,
-                                          const Array<const Vector *> &alfun,
-                                          const Array<const Vector *> &pelfun,
-                                          const Array<Vector *> &elvec)
-    {
-        int dof_u0 = el[0]->GetDof();
-        int dof_r0 = pel[0]->GetDof();
-        int dof_r1 = pel[1]->GetDof();
-
-        int dim = el[0]->GetDim();
-
-        Vector& e0 = *(elvec[0]);
-        Vector& e1 = *(elvec[1]);
-
-        e0.SetSize(dof_r0);
-        e0=0.0;
-        e1.SetSize(dof_r1);
-        e1=0.0;
-
-        int spaceDim = Tr.GetSpaceDim();
-        if (dim != spaceDim)
-        {
-            mfem::mfem_error(" PrmBlockLSFEMDiffusion::AssembleElementVector"
-                             " is not defined on manifold meshes");
-        }
-
-        //shape functions
-        Vector shu0(dof_u0);
-        Vector shr0(dof_r0);
-        Vector shr1(dof_r1);
-        DenseMatrix dsu0(dof_u0,dim);
-        DenseMatrix B(dof_u0, 4);
-        B=0.0;
-
-        double w;
-
-        Vector param(2); param=0.0;
-        Vector uu(4); uu=0.0;
-        Vector aa(4); aa=0.0;
-        Vector rr;
-        Vector lvec0; lvec0.SetSize(dof_r0);
-        Vector lvec1; lvec1.SetSize(dof_r1);
-
-        const IntegrationRule *ir = nullptr;
-        if(ir==nullptr){
-            int order= 2 * el[0]->GetOrder() + Tr.OrderGrad(el[0])
-                                +pel[0]->GetOrder()+pel[1]->GetOrder();
-            ir=&IntRules.Get(Tr.GetGeometryType(),order);
-        }
-
-        for (int i = 0; i < ir->GetNPoints(); i++)
-        {
-            const IntegrationPoint &ip = ir->IntPoint(i);
-            Tr.SetIntPoint(&ip);
-            w=Tr.Weight();
-            w = ip.weight * w;
-
-            el[0]->CalcPhysDShape(Tr,dsu0);
-            el[0]->CalcPhysShape(Tr,shu0);
-            pel[0]->CalcPhysShape(Tr,shr0);
-            pel[1]->CalcPhysShape(Tr,shr1);
-
-            param[0]=shr0*(*pelfun[0]);
-            param[1]=shr1*(*pelfun[1]);
-
-            //set the matrix B
-            for(int jj=0;jj<dim;jj++)
-            {
-                B.SetCol(jj,dsu0.GetColumn(jj));
-            }
-            B.SetCol(3,shu0);
-            B.MultTranspose(*elfun[0],uu);
-            B.MultTranspose(*alfun[0],aa);
-
-            qfunc->AQResidual(Tr, ip, param, uu, aa, rr);
-
-            lvec0=shr0;
-            lvec0*=rr[0];
-            lvec1=shr1;
-            lvec1*=rr[1];
-
-            e0.Add(w,lvec0);
-            e1.Add(w,lvec1);
-        }
-
-    }
-
-    virtual void AssemblePrmFaceVector(const Array<const FiniteElement *> &el1,
-                                       const Array<const FiniteElement *> &el2,
-                                       const Array<const FiniteElement *> &pel1,
-                                       const Array<const FiniteElement *> &pel2,
-                                       FaceElementTransformations &Tr,
-                                       const Array<const Vector *> &elfun,
-                                       const Array<const Vector *> &alfun,
-                                       const Array<const Vector *> &pelfun,
-                                       const Array<Vector *> &elvect)
-    {
-
-    }
-
-
-
-private:
-    LinDiffQFunc* qfunc;
-};
-
-}
-
-int main(int argc, char *argv[])
-{
-    const char *mesh_file = "../../data/beam-tet.mesh";
-    int ser_ref_levels = 1;
-    int order = 2;
-    bool visualization = true;
-    double newton_rel_tol = 1e-4;
-    double newton_abs_tol = 1e-6;
-    int newton_iter = 10;
-    int print_level = 0;
-
-    mfem::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(&order,
-                      "-o",
-                      "--order",
-                      "Order (degree) of the finite elements.");
-    args.AddOption(&visualization,
-                      "-vis",
-                      "--visualization",
-                      "-no-vis",
-                      "--no-visualization",
-                      "Enable or disable GLVis visualization.");
-    args.AddOption(&newton_rel_tol,
-                      "-rel",
-                      "--relative-tolerance",
-                      "Relative tolerance for the Newton solve.");
-    args.AddOption(&newton_abs_tol,
-                      "-abs",
-                      "--absolute-tolerance",
-                      "Absolute tolerance for the Newton solve.");
-    args.AddOption(&newton_iter,
-                      "-it",
-                      "--newton-iterations",
-                      "Maximum iterations for the Newton solve.");
-
-    args.Parse();
-    if (!args.Good())
-    {
-        args.PrintUsage(std::cout);
-        return 1;
-    }
-    args.PrintOptions(std::cout);
-    // 3. Read the (serial) mesh from the given mesh file on all processors.  We
-    //    can handle triangular, quadrilateral, tetrahedral and hexahedral meshes
-    //    with the same code.
-    mfem::Mesh *mesh = new mfem::Mesh(mesh_file, 1, 1);
-    int dim = mesh->Dimension();
-
-    // 4. Refine the mesh in serial to increase the resolution. In this example
-    //    we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
-    //    a command-line parameter.
-    for (int lev = 0; lev < ser_ref_levels; lev++)
-    {
-        mesh->UniformRefinement();
-    }
-
-    /// Define the q-function
-    mfem::ConstantCoefficient* dc=new mfem::ConstantCoefficient(1.0);
-    mfem::ConstantCoefficient* lc=new mfem::ConstantCoefficient(1.0);
-    mfem::LinDiffQFunc* qfun=new mfem::LinDiffQFunc(*dc,*lc,1.0,1.0,1.0);
-
-    mfem::H1_FECollection fec00(order, dim);
-    mfem::L2_FECollection fec01(order, dim);
-    mfem::FiniteElementSpace* bfes00=new mfem::FiniteElementSpace(mesh,&fec00,1,mfem::Ordering::byVDIM);
-    mfem::FiniteElementSpace* pfes00=new mfem::FiniteElementSpace(mesh,&fec00,1,mfem::Ordering::byVDIM);
-    mfem::FiniteElementSpace* pfes01=new mfem::FiniteElementSpace(mesh,&fec01,1,mfem::Ordering::byVDIM);
-
-    /// Define parametric nonlinear form
-    mfem::Array<mfem::FiniteElementSpace*> bfes;
-    mfem::Array<mfem::FiniteElementSpace*> pfes;
-
-    bfes.Append(bfes00);
-    pfes.Append(pfes00);
-    pfes.Append(pfes01);
-
-    mfem::PrmBlockNonlinearForm* nf=new mfem::PrmBlockNonlinearForm(bfes,pfes);
-    nf->AddDomainIntegrator(new mfem::PrmBlockLSFEMDiffusion(*qfun));
-
-    /// Define the grid functions
-    mfem::GridFunction* bgf00=new mfem::GridFunction(bfes00);
-    mfem::GridFunction* pgf00=new mfem::GridFunction(pfes00);
-    mfem::GridFunction* pgf01=new mfem::GridFunction(pfes01);
-    mfem::GridFunction* ggf00=new mfem::GridFunction(pfes00);
-    mfem::GridFunction* ggf01=new mfem::GridFunction(pfes01);
-
-    *bgf00=0.0;
-    *pgf00=1.0;
-    *pgf01=1.0;
-
-    mfem::BlockVector solbv; solbv.Update(nf->GetBlockTrueOffsets()); solbv=0.0;
-    mfem::BlockVector resbv; resbv.Update(nf->GetBlockTrueOffsets()); resbv=0.0;
-    mfem::BlockVector adjbv; adjbv.Update(nf->GetBlockTrueOffsets()); adjbv=0.0;
-    mfem::BlockVector prmbv; prmbv.Update(nf->PrmGetBlockTrueOffsets()); prmbv=1.0;
-    mfem::BlockVector grdbv; grdbv.Update(nf->PrmGetBlockTrueOffsets()); grdbv=0.0;
-
-    bgf00->SetFromTrueDofs(solbv.GetBlock(0));
-    pgf00->SetFromTrueDofs(prmbv.GetBlock(0));
-    pgf01->SetFromTrueDofs(prmbv.GetBlock(1));
-
-    nf->SetPrmFields(prmbv);
-    double energy = nf->GetEnergy(solbv);
-
-    nf->Mult(solbv,resbv);
-    std::cout<<"Norm res="<<resbv.Norml2()<<std::endl;
-
-    //mfem::Operator& K=nf->GetGradient(solbv);
-    std::cout<<"energy ="<< energy<<std::endl;
-
-    nf->SetStateFields(solbv);
-    nf->SetAdjointFields(adjbv);
-    nf->PrmMult(prmbv,grdbv);
-
-
-    //set the BC for the physics
-    mfem::Array<mfem::Array<int> *> ess_bdr;
-    mfem::Array<mfem::Vector*>      ess_rhs;
-    ess_bdr.Append(new mfem::Array<int>(mesh->bdr_attributes.Max()));
-    ess_rhs.Append(nullptr);
-    (*ess_bdr[0]) = 1;
-    nf->SetEssentialBC(ess_bdr,ess_rhs);
-
-    //define the solvers
-    mfem::UMFPackSolver* umfsolv=new mfem::UMFPackSolver();
-
-    mfem::GMRESSolver *gmres;
-    gmres = new mfem::GMRESSolver();
-    gmres->SetAbsTol(newton_abs_tol/10);
-    gmres->SetRelTol(newton_rel_tol/10);
-    gmres->SetMaxIter(100);
-    gmres->SetPrintLevel(print_level);
-    //gmres->SetPreconditioner(*prec);
-
-
-    mfem::NewtonSolver *ns;
-    ns = new mfem::NewtonSolver();
-    ns->iterative_mode = true;
-    ns->SetSolver(*gmres);
-    ns->SetOperator(*nf);
-    ns->SetPrintLevel(print_level);
-    ns->SetRelTol(newton_rel_tol);
-    ns->SetAbsTol(newton_abs_tol);
-    ns->SetMaxIter(newton_iter);
-
-    mfem::Vector b; //RHS is zero
-    solbv=0.0;
-    ns->Mult(b, solbv);
-
-
-    nf->SetStateFields(solbv);
-    nf->SetAdjointFields(solbv);
-    nf->PrmMult(prmbv,grdbv);
-
-
-
-
-    mfem::ParaViewDataCollection *dacol = new mfem::ParaViewDataCollection("Example91",
-                                                                      mesh);
-
-    ggf00->SetFromTrueDofs(grdbv.GetBlock(0));
-    ggf01->SetFromTrueDofs(grdbv.GetBlock(1));
-    pgf00->SetFromTrueDofs(solbv.GetBlock(0));
-
-    dacol->SetLevelsOfDetail(order);
-    dacol->RegisterField("sol", pgf00);
-    dacol->RegisterField("grad00", ggf00);
-    dacol->RegisterField("grad01", ggf01);
-
-    dacol->SetTime(1.0);
-    dacol->SetCycle(1);
-    dacol->Save();
-
-    delete dacol;
-
-    delete ns;
-    delete umfsolv;
-    delete gmres;
-    delete ess_bdr[0];
-
-    delete bgf00;
-    delete pgf00;
-    delete pgf01;
-    delete ggf00;
-    delete ggf01;
-
-
-    delete nf;
-    delete pfes01;
-    delete pfes00;
-    delete bfes00;
-
-    delete qfun;
-    delete lc;
-    delete dc;
-
-    delete mesh;
-
-}

examples/ex9p.cpp

@@ -64,66 +64,6 @@ double inflow_function(const Vector &x);
 // Mesh bounding box
 Vector bb_min, bb_max;
 
-// Type of preconditioner for implicit time integrator
-enum class PrecType : int
-{
-   ILU = 0,
-   AIR = 1
-};
-
-#if MFEM_HYPRE_VERSION >= 21800
-// Algebraic multigrid preconditioner for advective problems based on
-// approximate ideal restriction (AIR). Most effective when matrix is
-// first scaled by DG block inverse, and AIR applied to scaled matrix.
-// See https://doi.org/10.1137/17M1144350.
-class AIR_prec : public Solver
-{
-private:
-   const HypreParMatrix *A;
-   // Copy of A scaled by block-diagonal inverse
-   HypreParMatrix A_s;
-
-   HypreBoomerAMG *AIR_solver;
-   int blocksize;
-
-public:
-   AIR_prec(int blocksize_) : AIR_solver(NULL), blocksize(blocksize_) { }
-
-   void SetOperator(const Operator &op)
-   {
-      width = op.Width();
-      height = op.Height();
-
-      A = dynamic_cast<const HypreParMatrix *>(&op);
-      MFEM_VERIFY(A != NULL, "AIR_prec requires a HypreParMatrix.")
-
-      // Scale A by block-diagonal inverse
-      BlockInverseScale(A, &A_s, NULL, NULL, blocksize,
-                        BlockInverseScaleJob::MATRIX_ONLY);
-      delete AIR_solver;
-      AIR_solver = new HypreBoomerAMG(A_s);
-      AIR_solver->SetAdvectiveOptions(1, "", "FA");
-      AIR_solver->SetPrintLevel(0);
-      AIR_solver->SetMaxLevels(50);
-   }
-
-   virtual void Mult(const Vector &x, Vector &y) const
-   {
-      // Scale the rhs by block inverse and solve system
-      HypreParVector z_s;
-      BlockInverseScale(A, NULL, &x, &z_s, blocksize,
-                        BlockInverseScaleJob::RHS_ONLY);
-      AIR_solver->Mult(z_s, y);
-   }
-
-   ~AIR_prec()
-   {
-      delete AIR_solver;
-   }
-};
-#endif
-
-
 class DG_Solver : public Solver
 {
 private:
@@ -131,37 +71,24 @@ private:
    SparseMatrix M_diag;
    HypreParMatrix *A;
    GMRESSolver linear_solver;
-   Solver *prec;
+   BlockILU prec;
    double dt;
 public:
-   DG_Solver(HypreParMatrix &M_, HypreParMatrix &K_, const FiniteElementSpace &fes,
-             PrecType prec_type)
+   DG_Solver(HypreParMatrix &M_, HypreParMatrix &K_, const FiniteElementSpace &fes)
       : M(M_),
         K(K_),
         A(NULL),
         linear_solver(M.GetComm()),
+        prec(fes.GetFE(0)->GetDof(),
+             BlockILU::Reordering::MINIMUM_DISCARDED_FILL),
         dt(-1.0)
    {
-      int block_size = fes.GetFE(0)->GetDof();
-      if (prec_type == PrecType::ILU)
-      {
-         prec = new BlockILU(block_size,
-                             BlockILU::Reordering::MINIMUM_DISCARDED_FILL);
-      }
-      else if (prec_type == PrecType::AIR)
-      {
-#if MFEM_HYPRE_VERSION >= 21800
-         prec = new AIR_prec(block_size);
-#else
-         MFEM_ABORT("Must have MFEM_HYPRE_VERSION >= 21800 to use AIR.\n");
-#endif
-      }
       linear_solver.iterative_mode = false;
       linear_solver.SetRelTol(1e-9);
       linear_solver.SetAbsTol(0.0);
       linear_solver.SetMaxIter(100);
       linear_solver.SetPrintLevel(0);
-      linear_solver.SetPreconditioner(*prec);
+      linear_solver.SetPreconditioner(prec);
 
       M.GetDiag(M_diag);
    }
@@ -194,12 +121,10 @@ public:
 
    ~DG_Solver()
    {
-      delete prec;
       delete A;
    }
 };
 
-
 /** A time-dependent operator for the right-hand side of the ODE. The DG weak
     form of du/dt = -v.grad(u) is M du/dt = K u + b, where M and K are the mass
     and advection matrices, and b describes the flow on the boundary. This can
@@ -217,8 +142,7 @@ private:
    mutable Vector z;
 
 public:
-   FE_Evolution(ParBilinearForm &_M, ParBilinearForm &_K, const Vector &_b,
-                PrecType prec_type);
+   FE_Evolution(ParBilinearForm &_M, ParBilinearForm &_K, const Vector &_b);
 
    virtual void Mult(const Vector &x, Vector &y) const;
    virtual void ImplicitSolve(const double dt, const Vector &x, Vector &k);
@@ -230,9 +154,10 @@ public:
 int main(int argc, char *argv[])
 {
    // 1. Initialize MPI.
-   MPI_Session mpi;
-   int num_procs = mpi.WorldSize();
-   int myid = mpi.WorldRank();
+   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.
    problem = 0;
@@ -253,11 +178,7 @@ int main(int argc, char *argv[])
    bool adios2 = false;
    bool binary = false;
    int vis_steps = 5;
-#if MFEM_HYPRE_VERSION >= 21800
-   PrecType prec_type = PrecType::AIR;
-#else
-   PrecType prec_type = PrecType::ILU;
-#endif
+
    int precision = 8;
    cout.precision(precision);
 
@@ -291,8 +212,6 @@ int main(int argc, char *argv[])
                   "Final time; start time is 0.");
    args.AddOption(&dt, "-dt", "--time-step",
                   "Time step.");
-   args.AddOption((int *)&prec_type, "-pt", "--prec-type", "Preconditioner for "
-                  "implicit solves. 0 for ILU, 1 for pAIR-AMG.");
    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                   "--no-visualization",
                   "Enable or disable GLVis visualization.");
@@ -313,19 +232,20 @@ int main(int argc, char *argv[])
    args.Parse();
    if (!args.Good())
    {
-      if (mpi.Root())
+      if (myid == 0)
       {
          args.PrintUsage(cout);
       }
+      MPI_Finalize();
       return 1;
    }
-   if (mpi.Root())
+   if (myid == 0)
    {
       args.PrintOptions(cout);
    }
 
    Device device(device_config);
-   if (mpi.Root()) { device.Print(); }
+   if (myid == 0) { device.Print(); }
 
    // 3. Read the serial mesh from the given mesh file on all processors. We can
    //    handle geometrically periodic meshes in this code.
@@ -352,11 +272,12 @@ int main(int argc, char *argv[])
       case 23: ode_solver = new SDIRK23Solver; break;
       case 24: ode_solver = new SDIRK34Solver; break;
       default:
-         if (mpi.Root())
+         if (myid == 0)
          {
             cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
          }
          delete mesh;
+         MPI_Finalize();
          return 3;
    }
 
@@ -390,7 +311,7 @@ int main(int argc, char *argv[])
    ParFiniteElementSpace *fes = new ParFiniteElementSpace(pmesh, &fec);
 
    HYPRE_Int global_vSize = fes->GlobalTrueVSize();
-   if (mpi.Root())
+   if (myid == 0)
    {
       cout << "Number of unknowns: " << global_vSize << endl;
    }
@@ -421,16 +342,15 @@ int main(int argc, char *argv[])
    }
 
    m->AddDomainIntegrator(new MassIntegrator);
-   constexpr double alpha = -1.0;
-   k->AddDomainIntegrator(new ConvectionIntegrator(velocity, alpha));
+   k->AddDomainIntegrator(new ConvectionIntegrator(velocity, -1.0));
    k->AddInteriorFaceIntegrator(
-      new NonconservativeDGTraceIntegrator(velocity, alpha));
+      new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
    k->AddBdrFaceIntegrator(
-      new NonconservativeDGTraceIntegrator(velocity, alpha));
+      new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
 
    ParLinearForm *b = new ParLinearForm(fes);
    b->AddBdrFaceIntegrator(
-      new BoundaryFlowIntegrator(inflow, velocity, alpha));
+      new BoundaryFlowIntegrator(inflow, velocity, -1.0, -0.5));
 
    int skip_zeros = 0;
    m->Assemble();
@@ -531,11 +451,11 @@ int main(int argc, char *argv[])
       sout.open(vishost, visport);
       if (!sout)
       {
-         if (mpi.Root())
+         if (myid == 0)
             cout << "Unable to connect to GLVis server at "
                  << vishost << ':' << visport << endl;
          visualization = false;
-         if (mpi.Root())
+         if (myid == 0)
          {
             cout << "GLVis visualization disabled.\n";
          }
@@ -547,7 +467,7 @@ int main(int argc, char *argv[])
          sout << "solution\n" << *pmesh << *u;
          sout << "pause\n";
          sout << flush;
-         if (mpi.Root())
+         if (myid == 0)
             cout << "GLVis visualization paused."
                  << " Press space (in the GLVis window) to resume it.\n";
       }
@@ -556,7 +476,7 @@ int main(int argc, char *argv[])
    // 10. Define the time-dependent evolution operator describing the ODE
    //     right-hand side, and perform time-integration (looping over the time
    //     iterations, ti, with a time-step dt).
-   FE_Evolution adv(*m, *k, *B, prec_type);
+   FE_Evolution adv(*m, *k, *B);
 
    double t = 0.0;
    adv.SetTime(t);
@@ -573,7 +493,7 @@ int main(int argc, char *argv[])
 
       if (done || ti % vis_steps == 0)
       {
-         if (mpi.Root())
+         if (myid == 0)
          {
             cout << "time step: " << ti << ", time: " << t << endl;
          }
@@ -644,14 +564,16 @@ int main(int argc, char *argv[])
 #endif
    delete dc;
 
+   MPI_Finalize();
    return 0;
 }
 
 
 // Implementation of class FE_Evolution
 FE_Evolution::FE_Evolution(ParBilinearForm &_M, ParBilinearForm &_K,
-                           const Vector &_b, PrecType prec_type)
-   : TimeDependentOperator(_M.Height()), b(_b),
+                           const Vector &_b)
+   : TimeDependentOperator(_M.Height()),
+     b(_b),
      M_solver(_M.ParFESpace()->GetComm()),
      z(_M.Height())
 {
@@ -676,7 +598,7 @@ FE_Evolution::FE_Evolution(ParBilinearForm &_M, ParBilinearForm &_K,
       HypreSmoother *hypre_prec = new HypreSmoother(M_mat, HypreSmoother::Jacobi);
       M_prec = hypre_prec;
 
-      dg_solver = new DG_Solver(M_mat, K_mat, *_M.FESpace(), prec_type);
+      dg_solver = new DG_Solver(M_mat, K_mat, *_M.FESpace());
    }
    else
    {
@@ -692,10 +614,6 @@ FE_Evolution::FE_Evolution(ParBilinearForm &_M, ParBilinearForm &_K,
    M_solver.SetPrintLevel(0);
 }
 
-// Solve the equation:
-//    u_t = M^{-1}(Ku + b),
-// by solving associated linear system
-//    (M - dt*K) d = K*u + b
 void FE_Evolution::ImplicitSolve(const double dt, const Vector &x, Vector &k)
 {
    K->Mult(x, z);

fem/CMakeLists.txt

@@ -66,7 +66,6 @@ set(SRCS
   tmop_tools.cpp
   gslib.cpp
   transfer.cpp
-  prmnonlinearform.cpp
   )
 
 set(HDRS
@@ -113,7 +112,6 @@ set(HDRS
   tmop_tools.hpp
   gslib.hpp
   transfer.hpp
-  prmnonlinearform.hpp
   )
 
 if (MFEM_USE_SIDRE)
@@ -138,8 +136,7 @@ if (MFEM_USE_MPI)
     pgridfunc.cpp
     plinearform.cpp
     pnonlinearform.cpp
-    prestriction.cpp
-    pprmnonlinearform.cpp)
+    prestriction.cpp)
   # If this list (HDRS -> HEADERS) is used for install, we probably want the
   # headers added all the time.
   list(APPEND HDRS
@@ -148,8 +145,7 @@ if (MFEM_USE_MPI)
     pgridfunc.hpp
     plinearform.hpp
     pnonlinearform.hpp
-    prestriction.hpp
-    pprmnonlinearform.hpp)
+    prestriction.hpp)
 endif()
 
 convert_filenames_to_full_paths(SRCS)

fem/bilinearform.hpp

@@ -854,9 +854,9 @@ public:
 
       This returns the same operator as FormRectangularLinearSystem(), but does
       without the transformations of the right-hand side. */
-   virtual void FormRectangularSystemMatrix(const Array<int> &trial_tdof_list,
-                                            const Array<int> &test_tdof_list,
-                                            OperatorHandle &A);
+   void FormRectangularSystemMatrix(const Array<int> &trial_tdof_list,
+                                    const Array<int> &test_tdof_list,
+                                    OperatorHandle &A);
 
    /** @brief Form the column-constrained linear system matrix A.
        See FormRectangularSystemMatrix() for details.
@@ -883,11 +883,10 @@ public:
        Return in @a A a *reference* to the system matrix that is column-constrained.
        The reference will be invalidated when SetOperatorType(), Update(), or the
        destructor is called. */
-   virtual void FormRectangularLinearSystem(const Array<int> &trial_tdof_list,
-                                            const Array<int> &test_tdof_list,
-                                            Vector &x, Vector &b,
-                                            OperatorHandle &A, Vector &X,
-                                            Vector &B);
+   void FormRectangularLinearSystem(const Array<int> &trial_tdof_list,
+                                    const Array<int> &test_tdof_list,
+                                    Vector &x, Vector &b,
+                                    OperatorHandle &A, Vector &X, Vector &B);
 
    /** @brief Form the linear system A X = B, corresponding to this bilinear
        form and the linear form @a b(.).

fem/bilininteg.cpp

@@ -229,159 +229,6 @@ void SumIntegrator::AssembleElementMatrix(
    }
 }
 
-void SumIntegrator::AssembleElementMatrix2(
-   const FiniteElement &el1, const FiniteElement &el2,
-   ElementTransformation &Trans, DenseMatrix &elmat)
-{
-   MFEM_ASSERT(integrators.Size() > 0, "empty SumIntegrator.");
-
-   integrators[0]->AssembleElementMatrix2(el1, el2, Trans, elmat);
-   for (int i = 1; i < integrators.Size(); i++)
-   {
-      integrators[i]->AssembleElementMatrix2(el1, el2, Trans, elem_mat);
-      elmat += elem_mat;
-   }
-}
-
-void SumIntegrator::AssembleFaceMatrix(
-   const FiniteElement &el1, const FiniteElement &el2,
-   FaceElementTransformations &Trans, DenseMatrix &elmat)
-{
-   MFEM_ASSERT(integrators.Size() > 0, "empty SumIntegrator.");
-
-   integrators[0]->AssembleFaceMatrix(el1, el2, Trans, elmat);
-   for (int i = 1; i < integrators.Size(); i++)
-   {
-      integrators[i]->AssembleFaceMatrix(el1, el2, Trans, elem_mat);
-      elmat += elem_mat;
-   }
-}
-
-void SumIntegrator::AssembleFaceMatrix(
-   const FiniteElement &tr_fe,
-   const FiniteElement &te_fe1, const FiniteElement &te_fe2,
-   FaceElementTransformations &Trans, DenseMatrix &elmat)
-{
-   MFEM_ASSERT(integrators.Size() > 0, "empty SumIntegrator.");
-
-   integrators[0]->AssembleFaceMatrix(tr_fe, te_fe1, te_fe2, Trans, elmat);
-   for (int i = 1; i < integrators.Size(); i++)
-   {
-      integrators[i]->AssembleFaceMatrix(tr_fe, te_fe1, te_fe2, Trans, elem_mat);
-      elmat += elem_mat;
-   }
-}
-
-void SumIntegrator::AssemblePA(const FiniteElementSpace& fes)
-{
-   for (int i = 0; i < integrators.Size(); i++)
-   {
-      integrators[i]->AssemblePA(fes);
-   }
-}
-
-void SumIntegrator::AssembleDiagonalPA(Vector &diag)
-{
-   for (int i = 0; i < integrators.Size(); i++)
-   {
-      integrators[i]->AssembleDiagonalPA(diag);
-   }
-}
-
-void SumIntegrator::AssemblePAInteriorFaces(const FiniteElementSpace &fes)
-{
-   for (int i = 0; i < integrators.Size(); i++)
-   {
-      integrators[i]->AssemblePAInteriorFaces(fes);
-   }
-}
-
-void SumIntegrator::AssemblePABoundaryFaces(const FiniteElementSpace &fes)
-{
-   for (int i = 0; i < integrators.Size(); i++)
-   {
-      integrators[i]->AssemblePABoundaryFaces(fes);
-   }
-}
-
-void SumIntegrator::AddMultPA(const Vector& x, Vector& y) const
-{
-   for (int i = 0; i < integrators.Size(); i++)
-   {
-      integrators[i]->AddMultPA(x, y);
-   }
-}
-
-void SumIntegrator::AddMultTransposePA(const Vector &x, Vector &y) const
-{
-   for (int i = 0; i < integrators.Size(); i++)
-   {
-      integrators[i]->AddMultTransposePA(x, y);
-   }
-}
-
-void SumIntegrator::AssembleMF(const FiniteElementSpace &fes)
-{
-   for (int i = 0; i < integrators.Size(); i++)
-   {
-      integrators[i]->AssembleMF(fes);
-   }
-}
-
-void SumIntegrator::AddMultMF(const Vector& x, Vector& y) const
-{
-   for (int i = 0; i < integrators.Size(); i++)
-   {
-      integrators[i]->AddMultTransposeMF(x, y);
-   }
-}
-
-void SumIntegrator::AddMultTransposeMF(const Vector &x, Vector &y) const
-{
-   for (int i = 0; i < integrators.Size(); i++)
-   {
-      integrators[i]->AddMultMF(x, y);
-   }
-}
-
-void SumIntegrator::AssembleDiagonalMF(Vector &diag)
-{
-   for (int i = 0; i < integrators.Size(); i++)
-   {
-      integrators[i]->AssembleDiagonalMF(diag);
-   }
-}
-
-void SumIntegrator::AssembleEA(const FiniteElementSpace &fes, Vector &emat,
-                               const bool add)
-{
-   for (int i = 0; i < integrators.Size(); i++)
-   {
-      integrators[i]->AssembleEA(fes, emat, add);
-   }
-}
-
-void SumIntegrator::AssembleEAInteriorFaces(const FiniteElementSpace &fes,
-                                            Vector &ea_data_int,
-                                            Vector &ea_data_ext,
-                                            const bool add)
-{
-   for (int i = 0; i < integrators.Size(); i++)
-   {
-      integrators[i]->AssembleEAInteriorFaces(fes,ea_data_int,ea_data_ext,add);
-   }
-}
-
-void SumIntegrator::AssembleEABoundaryFaces(const FiniteElementSpace &fes,
-                                            Vector &ea_data_bdr,
-                                            const bool add)
-{
-   for (int i = 0; i < integrators.Size(); i++)
-   {
-      integrators[i]->AssembleEABoundaryFaces(fes, ea_data_bdr, add);
-   }
-}
-
 SumIntegrator::~SumIntegrator()
 {
    if (own_integrators)

fem/bilininteg.hpp

@@ -355,7 +355,7 @@ class SumIntegrator : public BilinearFormIntegrator
 {
 private:
    int own_integrators;
-   mutable DenseMatrix elem_mat;
+   DenseMatrix elem_mat;
    Array<BilinearFormIntegrator*> integrators;
 
 public:
@@ -367,55 +367,6 @@ public:
    virtual void AssembleElementMatrix(const FiniteElement &el,
                                       ElementTransformation &Trans,
                                       DenseMatrix &elmat);
-   virtual void AssembleElementMatrix2(const FiniteElement &trial_fe,
-                                       const FiniteElement &test_fe,
-                                       ElementTransformation &Trans,
-                                       DenseMatrix &elmat);
-
-   using BilinearFormIntegrator::AssembleFaceMatrix;
-   virtual void AssembleFaceMatrix(const FiniteElement &el1,
-                                   const FiniteElement &el2,
-                                   FaceElementTransformations &Trans,
-                                   DenseMatrix &elmat);
-
-   virtual void AssembleFaceMatrix(const FiniteElement &trial_face_fe,
-                                   const FiniteElement &test_fe1,
-                                   const FiniteElement &test_fe2,
-                                   FaceElementTransformations &Trans,
-                                   DenseMatrix &elmat);
-
-   using BilinearFormIntegrator::AssemblePA;
-   virtual void AssemblePA(const FiniteElementSpace& fes);
-
-   virtual void AssembleDiagonalPA(Vector &diag);
-
-   virtual void AssemblePAInteriorFaces(const FiniteElementSpace &fes);
-
-   virtual void AssemblePABoundaryFaces(const FiniteElementSpace &fes);
-
-   virtual void AddMultTransposePA(const Vector &x, Vector &y) const;
-
-   virtual void AddMultPA(const Vector& x, Vector& y) const;
-
-   virtual void AssembleMF(const FiniteElementSpace &fes);
-
-   virtual void AddMultMF(const Vector &x, Vector &y) const;
-
-   virtual void AddMultTransposeMF(const Vector &x, Vector &y) const;
-
-   virtual void AssembleDiagonalMF(Vector &diag);
-
-   virtual void AssembleEA(const FiniteElementSpace &fes, Vector &emat,
-                           const bool add);
-
-   virtual void AssembleEAInteriorFaces(const FiniteElementSpace &fes,
-                                        Vector &ea_data_int,
-                                        Vector &ea_data_ext,
-                                        const bool add);
-
-   virtual void AssembleEABoundaryFaces(const FiniteElementSpace &fes,
-                                        Vector &ea_data_bdr,
-                                        const bool add);
 
    virtual ~SumIntegrator();
 };
@@ -2045,8 +1996,6 @@ public:
 
    virtual void AddMultPA(const Vector&, Vector&) const;
 
-   virtual void AddMultTransposePA(const Vector&, Vector&) const;
-
    static const IntegrationRule &GetRule(const FiniteElement &trial_fe,
                                          const FiniteElement &test_fe);
 };
@@ -2110,8 +2059,6 @@ public:
 
    virtual void AddMultPA(const Vector&, Vector&) const;
 
-   virtual void AddMultTransposePA(const Vector&, Vector&) const;
-
    static const IntegrationRule &GetRule(const FiniteElement &trial_fe,
                                          const FiniteElement &test_fe,
                                          ElementTransformation &Trans);
@@ -2173,17 +2120,6 @@ public:
                                          ElementTransformation &Trans);
 };
 
-// Alias for @ConvectionIntegrator.
-using NonconservativeConvectionIntegrator = ConvectionIntegrator;
-
-/// -alpha (u, q . grad v), negative transpose of ConvectionIntegrator
-class ConservativeConvectionIntegrator : public TransposeIntegrator
-{
-public:
-   ConservativeConvectionIntegrator(VectorCoefficient &q, double a = 1.0)
-      : TransposeIntegrator(new ConvectionIntegrator(q, -a)) { }
-};
-
 /// alpha (q . grad u, v) using the "group" FE discretization
 class GroupConvectionIntegrator : public BilinearFormIntegrator
 {
@@ -2748,15 +2684,15 @@ public:
     One use case for this integrator is to discretize the operator -u.grad(v)
     with a DG formulation. The resulting formulation uses the
     ConvectionIntegrator (with coefficient u, and parameter alpha = -1) and the
-    transpose of the DGTraceIntegrator (with coefficient u, and parameters alpha
-    = 1, beta = -1/2 to use the upwind face flux, see also
-    NonconservativeDGTraceIntegrator). This discretization and the handling of
-    the inflow and outflow boundaries is illustrated in Example 9/9p.
+    transpose of the DGTraceIntegrator (with coefficient u, and parameters
+    alpha = 1, beta = -1/2 to use the upwind face flux). This discretization and
+    the handling of the inflow and outflow boundaries is illustrated in Example
+    9/9p.
 
     Another use case for this integrator is to discretize the operator -div(u v)
     with a DG formulation. The resulting formulation is conservative and
-    consists of the ConservativeConvectionIntegrator (with coefficient u, and
-    parameter alpha = -1) plus the DGTraceIntegrator (with coefficient u, and
+    consists of the transpose of the ConvectionIntegrator (with coefficient u,
+    and parameter alpha = 1) plus the DGTraceIntegrator (with coefficient u, and
     parameters alpha = -1, beta = -1/2 to use the upwind face flux).
     */
 class DGTraceIntegrator : public BilinearFormIntegrator
@@ -2775,17 +2711,13 @@ private:
    Vector shape1, shape2;
 
 public:
-   /// Construct integrator with rho = 1, b = 0.5*a.
-   DGTraceIntegrator(VectorCoefficient &u_, double a)
-   { rho = NULL; u = &u_; alpha = a; beta = 0.5*a; }
-
    /// Construct integrator with rho = 1.
-   DGTraceIntegrator(VectorCoefficient &u_, double a, double b)
-   { rho = NULL; u = &u_; alpha = a; beta = b; }
+   DGTraceIntegrator(VectorCoefficient &_u, double a, double b)
+   { rho = NULL; u = &_u; alpha = a; beta = b; }
 
-   DGTraceIntegrator(Coefficient &_rho, VectorCoefficient &u_,
+   DGTraceIntegrator(Coefficient &_rho, VectorCoefficient &_u,
                      double a, double b)
-   { rho = &_rho; u = &u_; alpha = a; beta = b; }
+   { rho = &_rho; u = &_u; alpha = a; beta = b; }
 
    using BilinearFormIntegrator::AssembleFaceMatrix;
    virtual void AssembleFaceMatrix(const FiniteElement &el1,
@@ -2819,30 +2751,6 @@ private:
    void SetupPA(const FiniteElementSpace &fes, FaceType type);
 };
 
-// Alias for @a DGTraceIntegrator.
-using ConservativeDGTraceIntegrator = DGTraceIntegrator;
-
-/** Integrator that represents the face terms used for the non-conservative
-    DG discretization of the convection equation:
-    -alpha < rho_u (u.n) {v},[w] > + beta < rho_u |u.n| [v],[w] >.
-
-    This integrator can be used with together with ConvectionIntegrator to
-    implement an upwind DG discretization in non-conservative form, see ex9 and
-    ex9p. */
-class NonconservativeDGTraceIntegrator : public TransposeIntegrator
-{
-public:
-   NonconservativeDGTraceIntegrator(VectorCoefficient &u, double a)
-      : TransposeIntegrator(new DGTraceIntegrator(u, -a, 0.5*a)) { }
-
-   NonconservativeDGTraceIntegrator(VectorCoefficient &u, double a, double b)
-      : TransposeIntegrator(new DGTraceIntegrator(u, -a, b)) { }
-
-   NonconservativeDGTraceIntegrator(Coefficient &rho, VectorCoefficient &u,
-                                    double a, double b)
-      : TransposeIntegrator(new DGTraceIntegrator(rho, u, -a, b)) { }
-};
-
 /** Integrator for the DG form:
 
     - < {(Q grad(u)).n}, [v] > + sigma < [u], {(Q grad(v)).n} >

fem/bilininteg_diffusion_pa.cpp

@@ -1904,17 +1904,4 @@ void DiffusionIntegrator::AddMultPA(const Vector &x, Vector &y) const
    }
 }
 
-void DiffusionIntegrator::AddMultTransposePA(const Vector &x, Vector &y) const
-{
-   if (symmetric)
-   {
-      AddMultPA(x, y);
-   }
-   else
-   {
-      MFEM_ABORT("DiffusionIntegrator::AddMultTransposePA only implemented in "
-                 "the symmetric case.")
-   }
-}
-
 } // namespace mfem

fem/bilininteg_mass_pa.cpp

@@ -1226,10 +1226,4 @@ void MassIntegrator::AddMultPA(const Vector &x, Vector &y) const
    }
 }
 
-void MassIntegrator::AddMultTransposePA(const Vector &x, Vector &y) const
-{
-   // Mass integrator is symmetric
-   AddMultPA(x, y);
-}
-
 } // namespace mfem

fem/fe.cpp

@@ -7797,10 +7797,7 @@ NodalTensorFiniteElement::NodalTensorFiniteElement(const int dims,
                                                    const DofMapType dmtype)
    : NodalFiniteElement(dims, GetTensorProductGeometry(dims), Pow(p + 1, dims),
                         p, dims > 1 ? FunctionSpace::Qk : FunctionSpace::Pk),
-     TensorBasisElement(dims, p, VerifyNodal(btype), dmtype)
-{
-   lex_ordering = dof_map;
-}
+     TensorBasisElement(dims, p, VerifyNodal(btype), dmtype) { }
 
 
 PositiveTensorFiniteElement::PositiveTensorFiniteElement(
@@ -8484,33 +8481,23 @@ H1_TriangleElement::H1_TriangleElement(const int p, const int btype)
    Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1);
 #endif
 
-   int p2p3 = 2*p + 3;
-   auto idx = [p2p3](int i, int j) { return ((p2p3-j)*j)/2+i; };
-   lex_ordering.SetSize(dof);
-
    // vertices
-   lex_ordering[idx(0,0)] = 0;
    Nodes.IntPoint(0).Set2(cp[0], cp[0]);
-   lex_ordering[idx(p,0)] = 1;
    Nodes.IntPoint(1).Set2(cp[p], cp[0]);
-   lex_ordering[idx(0,p)] = 2;
    Nodes.IntPoint(2).Set2(cp[0], cp[p]);
 
    // edges
    int o = 3;
    for (int i = 1; i < p; i++)
    {
-      lex_ordering[idx(i,0)] = o;
       Nodes.IntPoint(o++).Set2(cp[i], cp[0]);
    }
    for (int i = 1; i < p; i++)
    {
-      lex_ordering[idx(p-i,i)] = o;
       Nodes.IntPoint(o++).Set2(cp[p-i], cp[i]);
    }
    for (int i = 1; i < p; i++)
    {
-      lex_ordering[idx(0,p-i)] = o;
       Nodes.IntPoint(o++).Set2(cp[0], cp[p-i]);
    }
 
@@ -8519,7 +8506,6 @@ H1_TriangleElement::H1_TriangleElement(const int p, const int btype)
       for (int i = 1; i + j < p; i++)
       {
          const double w = cp[i] + cp[j] + cp[p-i-j];
-         lex_ordering[idx(i,j)] = o;
          Nodes.IntPoint(o++).Set2(cp[i]/w, cp[j]/w);
       }
 
@@ -8653,56 +8639,36 @@ H1_TetrahedronElement::H1_TetrahedronElement(const int p, const int btype)
    Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
 #endif
 
-   auto tri = [](int k) { return (k*(k + 1))/2; };
-   auto tet = [](int k) { return (k*(k + 1)*(k + 2))/6; };
-   int ndof = tet(p+1);
-   auto idx = [tri, tet, p, ndof](int i, int j, int k)
-   {
-      return ndof - tet(p - k) - tri(p + 1 - k - j) + i;
-   };
-
-   lex_ordering.SetSize(dof);
-
    // vertices
-   lex_ordering[idx(0,0,0)] = 0;
    Nodes.IntPoint(0).Set3(cp[0], cp[0], cp[0]);
-   lex_ordering[idx(p,0,0)] = 1;
    Nodes.IntPoint(1).Set3(cp[p], cp[0], cp[0]);
-   lex_ordering[idx(0,p,0)] = 2;
    Nodes.IntPoint(2).Set3(cp[0], cp[p], cp[0]);
-   lex_ordering[idx(0,0,p)] = 3;
    Nodes.IntPoint(3).Set3(cp[0], cp[0], cp[p]);
 
    // edges (see Tetrahedron::edges in mesh/tetrahedron.cpp)
    int o = 4;
    for (int i = 1; i < p; i++)  // (0,1)
    {
-      lex_ordering[idx(i,0,0)] = o;
       Nodes.IntPoint(o++).Set3(cp[i], cp[0], cp[0]);
    }
    for (int i = 1; i < p; i++)  // (0,2)
    {
-      lex_ordering[idx(0,i,0)] = o;
       Nodes.IntPoint(o++).Set3(cp[0], cp[i], cp[0]);
    }
    for (int i = 1; i < p; i++)  // (0,3)
    {
-      lex_ordering[idx(0,0,i)] = o;
       Nodes.IntPoint(o++).Set3(cp[0], cp[0], cp[i]);
    }
    for (int i = 1; i < p; i++)  // (1,2)
    {
-      lex_ordering[idx(p-i,i,0)] = o;
       Nodes.IntPoint(o++).Set3(cp[p-i], cp[i], cp[0]);
    }
    for (int i = 1; i < p; i++)  // (1,3)
    {
-      lex_ordering[idx(p-i,0,i)] = o;
       Nodes.IntPoint(o++).Set3(cp[p-i], cp[0], cp[i]);
    }
    for (int i = 1; i < p; i++)  // (2,3)
    {
-      lex_ordering[idx(0,p-i,i)] = o;
       Nodes.IntPoint(o++).Set3(cp[0], cp[p-i], cp[i]);
    }
 
@@ -8710,28 +8676,24 @@ H1_TetrahedronElement::H1_TetrahedronElement(const int p, const int btype)
    for (int j = 1; j < p; j++)
       for (int i = 1; i + j < p; i++)  // (1,2,3)
       {
-         lex_ordering[idx(p-i-j,i,j)] = o;
          double w = cp[i] + cp[j] + cp[p-i-j];
          Nodes.IntPoint(o++).Set3(cp[p-i-j]/w, cp[i]/w, cp[j]/w);
       }
    for (int j = 1; j < p; j++)
       for (int i = 1; i + j < p; i++)  // (0,3,2)
       {
-         lex_ordering[idx(0,j,i)] = o;
          double w = cp[i] + cp[j] + cp[p-i-j];
          Nodes.IntPoint(o++).Set3(cp[0], cp[j]/w, cp[i]/w);
       }
    for (int j = 1; j < p; j++)
       for (int i = 1; i + j < p; i++)  // (0,1,3)
       {
-         lex_ordering[idx(i,0,j)] = o;
          double w = cp[i] + cp[j] + cp[p-i-j];
          Nodes.IntPoint(o++).Set3(cp[i]/w, cp[0], cp[j]/w);
       }
    for (int j = 1; j < p; j++)
       for (int i = 1; i + j < p; i++)  // (0,2,1)
       {
-         lex_ordering[idx(j,i,0)] = o;
          double w = cp[i] + cp[j] + cp[p-i-j];
          Nodes.IntPoint(o++).Set3(cp[j]/w, cp[i]/w, cp[0]);
       }
@@ -8741,7 +8703,6 @@ H1_TetrahedronElement::H1_TetrahedronElement(const int p, const int btype)
       for (int j = 1; j + k < p; j++)
          for (int i = 1; i + j + k < p; i++)
          {
-            lex_ordering[idx(i,j,k)] = o;
             double w = cp[i] + cp[j] + cp[k] + cp[p-i-j-k];
             Nodes.IntPoint(o++).Set3(cp[i]/w, cp[j]/w, cp[k]/w);
          }
@@ -9286,22 +9247,7 @@ H1_WedgeElement::H1_WedgeElement(const int p,
    t_dof.SetSize(dof);
    s_dof.SetSize(dof);
 
-   int p2p3 = 2*p + 3, ntri = ((p + 1)*(p + 2))/2;
-   auto idx = [p2p3,ntri](int i, int j, int k)
-   {
-      return k*ntri + ((p2p3-j)*j)/2+i;
-   };
-
-   lex_ordering.SetSize(dof);
-   int o = 0;
-
    // Nodal DoFs
-   lex_ordering[idx(0,0,0)] = o++;
-   lex_ordering[idx(p,0,0)] = o++;
-   lex_ordering[idx(0,p,0)] = o++;
-   lex_ordering[idx(0,0,p)] = o++;
-   lex_ordering[idx(p,0,p)] = o++;
-   lex_ordering[idx(0,p,p)] = o++;
    t_dof[0] = 0; s_dof[0] = 0;
    t_dof[1] = 1; s_dof[1] = 0;
    t_dof[2] = 2; s_dof[2] = 0;
@@ -9310,19 +9256,9 @@ H1_WedgeElement::H1_WedgeElement(const int p,
    t_dof[5] = 2; s_dof[5] = 1;
 
    // Edge DoFs
-   int k = 0;
    int ne = p-1;
    for (int i=1; i<p; i++)
    {
-      lex_ordering[idx(i,0,0)] = o + 0*ne + k;
-      lex_ordering[idx(p-i,i,0)] = o + 1*ne + k;
-      lex_ordering[idx(0,p-i,0)] = o + 2*ne + k;
-      lex_ordering[idx(i,0,p)] = o + 3*ne + k;
-      lex_ordering[idx(p-i,i,p)] = o + 4*ne + k;
-      lex_ordering[idx(0,p-i,p)] = o + 5*ne + k;
-      lex_ordering[idx(0,0,i)] = o + 6*ne + k;
-      lex_ordering[idx(p,0,i)] = o + 7*ne + k;
-      lex_ordering[idx(0,p,i)] = o + 8*ne + k;
       t_dof[5 + 0 * ne + i] = 2 + 0 * ne + i; s_dof[5 + 0 * ne + i] = 0;
       t_dof[5 + 1 * ne + i] = 2 + 1 * ne + i; s_dof[5 + 1 * ne + i] = 0;
       t_dof[5 + 2 * ne + i] = 2 + 2 * ne + i; s_dof[5 + 2 * ne + i] = 0;
@@ -9332,26 +9268,21 @@ H1_WedgeElement::H1_WedgeElement(const int p,
       t_dof[5 + 6 * ne + i] = 0;              s_dof[5 + 6 * ne + i] = i + 1;
       t_dof[5 + 7 * ne + i] = 1;              s_dof[5 + 7 * ne + i] = i + 1;
       t_dof[5 + 8 * ne + i] = 2;              s_dof[5 + 8 * ne + i] = i + 1;
-      ++k;
    }
-   o += 9*ne;
 
    // Triangular Face DoFs
-   k=0;
+   int k=0;
    int nt = (p-1)*(p-2)/2;
    for (int j=1; j<p; j++)
    {
       for (int i=1; i<p-j; i++)
       {
          int l = j - p + (((2 * p - 1) - i) * i) / 2;
-         lex_ordering[idx(i,j,0)] = o+l;
-         lex_ordering[idx(i,j,p)] = o+nt+k;
          t_dof[6 + 9 * ne + k]      = 3 * p + l; s_dof[6 + 9 * ne + k]      = 0;
          t_dof[6 + 9 * ne + nt + k] = 3 * p + k; s_dof[6 + 9 * ne + nt + k] = 1;
          k++;
       }
    }
-   o += 2*nt;
 
    // Quadrilateral Face DoFs
    k=0;
@@ -9360,10 +9291,6 @@ H1_WedgeElement::H1_WedgeElement(const int p,
    {
       for (int i=1; i<p; i++)
       {
-         lex_ordering[idx(i,0,j)] = o+k;
-         lex_ordering[idx(p-i,i,j)] = o+nq+k;
-         lex_ordering[idx(0,p-i,j)] = o+2*nq+k;
-
          t_dof[6 + 9 * ne + 2 * nt + 0 * nq + k] = 2 + 0 * ne + i;
          t_dof[6 + 9 * ne + 2 * nt + 1 * nq + k] = 2 + 1 * ne + i;
          t_dof[6 + 9 * ne + 2 * nt + 2 * nq + k] = 2 + 2 * ne + i;
@@ -9375,7 +9302,6 @@ H1_WedgeElement::H1_WedgeElement(const int p,
          k++;
       }
    }
-   o += 3*nq;
 
    // Interior DoFs
    int m=0;
@@ -9384,9 +9310,8 @@ H1_WedgeElement::H1_WedgeElement(const int p,
       int l=0;
       for (int j=1; j<p; j++)
       {
-         for (int i=1; i+j<p; i++)
+         for (int i=1; i<j; i++)
          {
-            lex_ordering[idx(i,j,k)] = o++;
             t_dof[6 + 9 * ne + 2 * nt + 3 * nq + m] = 3 * p + l;
             s_dof[6 + 9 * ne + 2 * nt + 3 * nq + m] = 1 + k;
             l++; m++;

fem/fe.hpp

@@ -698,7 +698,6 @@ public:
 class NodalFiniteElement : public ScalarFiniteElement
 {
 protected:
-   Array<int> lex_ordering;
    void ProjectCurl_2D(const FiniteElement &fe,
                        ElementTransformation &Trans,
                        DenseMatrix &curl) const;
@@ -747,29 +746,6 @@ public:
    virtual void ProjectDiv(const FiniteElement &fe,
                            ElementTransformation &Trans,
                            DenseMatrix &div) const;
-
-   /** @brief Get an Array<int> that maps lexicographically ordered indices to
-       the indices of the respective nodes/dofs/basis functions. Lexicographic
-       ordering of nodes is defined in terms of reference-space coordinates
-       (x,y,z). Lexicographically ordered nodes are listed first in order of
-       increasing x-coordinate, and then in order of increasing y-coordinate,
-       and finally in order of increasing z-coordinate.
-
-       For example, the six nodes of a quadratic triangle are lexicographically
-       ordered as follows:
-
-       5
-       |\
-       3 4
-       |  \
-       0-1-2
-
-       The resulting array may be empty if the DOFs are already ordered
-       lexicographically, or if the finite element does not support creating
-       this permutation. The array returned is the same as the array given by
-       TensorBasisElement::GetDofMap, but it is also available for non-tensor
-       elements. */
-   const Array<int> &GetLexicographicOrdering() const { return lex_ordering; }
 };
 
 /** @brief Class for finite elements utilizing the

fem/fe_coll.cpp

@@ -1854,16 +1854,19 @@ const int *H1_FECollection::GetDofMap(Geometry::Type GeomType) const
 {
    const int *dof_map = NULL;
    const FiniteElement *fe = H1_Elements[GeomType];
-   const NodalFiniteElement *nodal_fe =
-      dynamic_cast<const NodalFiniteElement*>(fe);
-   if (nodal_fe)
+   switch (GeomType)
    {
-      dof_map = nodal_fe->GetLexicographicOrdering().GetData();
-   }
-   else
-   {
-      MFEM_ABORT("Geometry type " << Geometry::Name[GeomType] << " is not "
-                 "implemented");
+      case Geometry::SEGMENT:
+      case Geometry::SQUARE:
+      case Geometry::CUBE:
+         dof_map = dynamic_cast<const TensorBasisElement *>(fe)
+                   ->GetDofMap().GetData();
+         break;
+      default:
+         MFEM_ABORT("Geometry type " << Geometry::Name[GeomType] << " is not "
+                    "implemented");
+         // The "Cartesian" ordering for other geometries is defined by the
+         // class GeometryRefiner.
    }
    return dof_map;
 }

fem/fem.hpp

@@ -41,7 +41,6 @@
 #include "transfer.hpp"
 #include "fespacehierarchy.hpp"
 #include "multigrid.hpp"
-#include "prmnonlinearform.hpp"
 
 #ifdef MFEM_USE_MPI
 #include "pfespace.hpp"
@@ -49,7 +48,6 @@
 #include "plinearform.hpp"
 #include "pbilinearform.hpp"
 #include "pnonlinearform.hpp"
-#include "pprmnonlinearform.hpp"
 #endif
 
 #ifdef MFEM_USE_SIDRE

fem/geom.cpp

@@ -1165,31 +1165,24 @@ RefinedGeometry * GeometryRefiner::Refine(Geometry::Type Geom,
          Array<int> vi((n+1)*(n+1)*(n+1));
          vi = -1;
          m = 0;
-
-         // vertices are given in lexicographic ordering on the reference
-         // element
-         for (int kk = 0; kk <= n; kk++)
-            for (int jj = 0; jj <= n-kk; jj++)
-               for (int ii = 0; ii <= n-jj-kk; ii++)
+         for (k = 0; k <= n; k++)
+            for (j = 0; j <= k; j++)
+               for (i = 0; i <= j; i++)
                {
                   IntegrationPoint &ip = RG->RefPts.IntPoint(m);
-                  double w = cp[ii] + cp[jj] + cp[kk] + cp[Times-ii-jj-kk];
-                  ip.x = cp[ii]/w;
-                  ip.y = cp[jj]/w;
-                  ip.z = cp[kk]/w;
-                  // (ii,jj,kk) are coordinates in the reference tetrahedron,
-                  // transform to coordinates (i,j,k) in the auxiliary
-                  // tetrahedron defined by (0,0,0), (0,0,1), (1,1,1), (0,1,1)
-                  int i = jj;
-                  int j = jj+kk;
-                  int k = ii+jj+kk;
+                  // map the coordinates to the reference tetrahedron
+                  // (0,0,0) -> (0,0,0)
+                  // (0,0,1) -> (1,0,0)
+                  // (1,1,1) -> (0,1,0)
+                  // (0,1,1) -> (0,0,1)
+                  double w = cp[k-j] + cp[i] + cp[j-i] + cp[Times-k];
+                  ip.x = cp[k-j]/w;
+                  ip.y = cp[i]/w;
+                  ip.z = cp[j-i]/w;
                   l = i + (j + k * (n+1)) * (n+1);
-                  // map from linear Cartesian hex index in the auxiliary tet
-                  // to lexicographic in the reference tet
                   vi[l] = m;
                   m++;
                }
-
          if (m != (n+3)*(n+2)*(n+1)/6)
          {
             mfem_error("GeometryRefiner::Refine() for TETRAHEDRON #1");
@@ -1276,9 +1269,18 @@ RefinedGeometry * GeometryRefiner::Refine(Geometry::Type Geom,
                for (i = 0; i <= n-j; i++, l++)
                {
                   IntegrationPoint &ip = RG->RefPts.IntPoint(l);
-                  ip.x = cp[i]/(cp[i] + cp[j] + cp[n-i-j]);
-                  ip.y = cp[j]/(cp[i] + cp[j] + cp[n-i-j]);
-                  ip.z = cp[k];
+                  if (type == 0)
+                  {
+                     ip.x = double(i) / n;
+                     ip.y = double(j) / n;
+                     ip.z = double(k) / n;
+                  }
+                  else
+                  {
+                     ip.x = cp[i]/(cp[i] + cp[j] + cp[n-i-j]);
+                     ip.y = cp[j]/(cp[i] + cp[j] + cp[n-i-j]);
+                     ip.z = cp[k];
+                  }
                   m++;
                }
          if (m != (n+1)*(n+1)*(n+2)/2)

fem/lininteg.hpp

@@ -406,10 +406,6 @@ private:
    Vector shape;
 
 public:
-   BoundaryFlowIntegrator(Coefficient &_f, VectorCoefficient &_u,
-                          double a)
-   { f = &_f; u = &_u; alpha = a; beta = 0.5*a; }
-
    BoundaryFlowIntegrator(Coefficient &_f, VectorCoefficient &_u,
                           double a, double b)
    { f = &_f; u = &_u; alpha = a; beta = b; }

fem/nonlininteg.cpp

@@ -128,92 +128,6 @@ double BlockNonlinearFormIntegrator::GetElementEnergy(
    return 0.0;
 }
 
-double PrmBlockNonlinearFormIntegrator::GetElementEnergy(const Array<const FiniteElement *> &el,
-                                                         const Array<const FiniteElement *> &pel,
-                                                         ElementTransformation &Tr,
-                                                         const Array<const Vector *> &elfun,
-                                                         const Array<const Vector *> &pelfun)
-{
-    mfem_error("PrmBlockNonlinearFormIntegrator::GetElementEnergy"
-               " is not overloaded!");
-    return 0.0;
-}
-
-void PrmBlockNonlinearFormIntegrator::AssembleFaceGrad(const Array<const FiniteElement *> &el1,
-                                                       const Array<const FiniteElement *> &el2,
-                                                       const Array<const FiniteElement *> &pel1,
-                                                       const Array<const FiniteElement *> &pel2,
-                                                       FaceElementTransformations &Tr,
-                                                       const Array<const Vector *> &elfun,
-                                                       const Array<const Vector *> &pelfun,
-                                                       const Array2D<DenseMatrix *> &elmats)
-{
-    mfem_error("PrmBlockNonlinearFormIntegrator::AssembleFaceGrad"
-               " is not overloaded!");
-}
-
-void PrmBlockNonlinearFormIntegrator::AssembleElementGrad(const Array<const FiniteElement *> &el,
-                                                          const Array<const FiniteElement *> &pel,
-                                                          ElementTransformation &Tr,
-                                                          const Array<const Vector *> &elfun,
-                                                          const Array<const Vector *> &pelfun,
-                                                          const Array2D<DenseMatrix *> &elmats)
-{
-    mfem_error("PrmBlockNonlinearFormIntegrator::AssembleElementGrad"
-               " is not overloaded!");
-}
-
-void PrmBlockNonlinearFormIntegrator::AssembleElementVector(const Array<const FiniteElement *> &el,
-                                                            const Array<const FiniteElement *> &pel,
-                                                            ElementTransformation &Tr,
-                                                            const Array<const Vector *> &elfun,
-                                                            const Array<const Vector *> &pelfun,
-                                                            const Array<Vector *> &elvec)
-{
-    mfem_error("PrmBlockNonlinearFormIntegrator::AssembleElementVector"
-               " is not overloaded!");
-}
-
-void PrmBlockNonlinearFormIntegrator::AssembleFaceVector(const Array<const FiniteElement *> &el1,
-                                                         const Array<const FiniteElement *> &el2,
-                                                         const Array<const FiniteElement *> &pel1,
-                                                         const Array<const FiniteElement *> &pel2,
-                                                         FaceElementTransformations &Tr,
-                                                         const Array<const Vector *> &elfun,
-                                                         const Array<const Vector *> &pelfun,
-                                                         const Array<Vector *> &elvect)
-{
-    mfem_error("PrmBlockNonlinearFormIntegrator::AssembleFaceVector"
-               " is not overloaded!");
-}
-
-
-void PrmBlockNonlinearFormIntegrator::AssemblePrmElementVector(const Array<const FiniteElement *> &el,
-                                                               const Array<const FiniteElement *> &pel,
-                                                               ElementTransformation &Tr,
-                                                               const Array<const Vector *> &elfun,
-                                                               const Array<const Vector *> &alfun,
-                                                               const Array<const Vector *> &pelfun,
-                                                               const Array<Vector *> &elvec)
-{
-    mfem_error("PrmBlockNonlinearFormIntegrator::AssemblePrmElementVector"
-               " is not overloaded!");
-}
-
-void PrmBlockNonlinearFormIntegrator::AssemblePrmFaceVector(const Array<const FiniteElement *> &el1,
-                                                         const Array<const FiniteElement *> &el2,
-                                                         const Array<const FiniteElement *> &pel1,
-                                                         const Array<const FiniteElement *> &pel2,
-                                                         FaceElementTransformations &Tr,
-                                                         const Array<const Vector *> &elfun,
-                                                         const Array<const Vector *> &alfun,
-                                                         const Array<const Vector *> &pelfun,
-                                                         const Array<Vector *> &elvect)
-{
-    mfem_error("PrmBlockNonlinearFormIntegrator::AssemblePrmFaceVector"
-               " is not overloaded!");
-}
-
 
 double InverseHarmonicModel::EvalW(const DenseMatrix &J) const
 {

fem/nonlininteg.hpp

@@ -130,80 +130,6 @@ public:
 };
 
 
-/** The abstract base class PrmBlockNonlinearFormIntegrator is
-    a generalization of the BlockNonlinearFormIntegrator class suitable
-    for block state and parameter vectors. */
-class PrmBlockNonlinearFormIntegrator
-{
-public:
-   /// Compute the local energy
-   virtual double GetElementEnergy(const Array<const FiniteElement *>&el,
-                                   const Array<const FiniteElement *>&pel,
-                                   ElementTransformation &Tr,
-                                   const Array<const Vector *>&elfun,
-                                   const Array<const Vector *>&pelfun);
-
-   /// Perform the local action of the BlockNonlinearFormIntegrator
-   virtual void AssembleElementVector(const Array<const FiniteElement *> &el,
-                                      const Array<const FiniteElement *>&pel,
-                                      ElementTransformation &Tr,
-                                      const Array<const Vector *> &elfun,
-                                      const Array<const Vector *>&pelfun,
-                                      const Array<Vector *> &elvec);
-
-   virtual void AssembleFaceVector(const Array<const FiniteElement *> &el1,
-                                   const Array<const FiniteElement *> &el2,
-                                   const Array<const FiniteElement *> &pel1,
-                                   const Array<const FiniteElement *> &pel2,
-                                   FaceElementTransformations &Tr,
-                                   const Array<const Vector *> &elfun,
-                                   const Array<const Vector *>&pelfun,
-                                   const Array<Vector *> &elvect);
-
-    /// Perform the local action on the parameters of the BlockNonlinearFormIntegrator
-    virtual void AssemblePrmElementVector(const Array<const FiniteElement *> &el,
-                                       const Array<const FiniteElement *>&pel,
-                                       ElementTransformation &Tr,
-                                       const Array<const Vector *> &elfun,
-                                       const Array<const Vector *> &alfun,
-                                       const Array<const Vector *>&pelfun,
-                                       const Array<Vector *> &elvec);
-
-    virtual void AssemblePrmFaceVector(const Array<const FiniteElement *> &el1,
-                                    const Array<const FiniteElement *> &el2,
-                                    const Array<const FiniteElement *> &pel1,
-                                    const Array<const FiniteElement *> &pel2,
-                                    FaceElementTransformations &Tr,
-                                    const Array<const Vector *> &elfun,
-                                    const Array<const Vector *> &alfun,
-                                    const Array<const Vector *>&pelfun,
-                                    const Array<Vector *> &elvect);
-
-
-
-   /// Assemble the local gradient matrix
-   virtual void AssembleElementGrad(const Array<const FiniteElement*> &el,
-                                    const Array<const FiniteElement *>&pel,
-                                    ElementTransformation &Tr,
-                                    const Array<const Vector *> &elfun,
-                                    const Array<const Vector *>&pelfun,
-                                    const Array2D<DenseMatrix *> &elmats);
-
-   virtual void AssembleFaceGrad(const Array<const FiniteElement *>&el1,
-                                 const Array<const FiniteElement *>&el2,
-                                 const Array<const FiniteElement *> &pel1,
-                                 const Array<const FiniteElement *> &pel2,
-                                 FaceElementTransformations &Tr,
-                                 const Array<const Vector *> &elfun,
-                                 const Array<const Vector *>&pelfun,
-                                 const Array2D<DenseMatrix *> &elmats);
-
-
-   virtual ~PrmBlockNonlinearFormIntegrator() { }
-};
-
-
-
 /// Abstract class for hyperelastic models
 class HyperelasticModel
 {

fem/pbilinearform.hpp

@@ -249,9 +249,6 @@ public:
        @a A. */
    void ParallelAssemble(OperatorHandle &A);
 
-   using MixedBilinearForm::FormRectangularSystemMatrix;
-   using MixedBilinearForm::FormRectangularLinearSystem;
-
    /** @brief Return in @a A a parallel (on truedofs) version of this operator.
 
        This returns the same operator as FormRectangularLinearSystem(), but does
@@ -316,8 +313,6 @@ public:
        domain and range parallel finite element spaces */
    void GetParBlocks(Array2D<HypreParMatrix *> &blocks) const;
 
-   using MixedBilinearForm::FormRectangularSystemMatrix;
-
    /** @brief Return in @a A a parallel (on truedofs) version of this operator. */
    virtual void FormRectangularSystemMatrix(OperatorHandle &A);
 

fem/pprmnonlinearform.cpp

@@ -1,363 +0,0 @@
-// Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
-// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
-// LICENSE and NOTICE for details. LLNL-CODE-806117.
-//
-// This file is part of the MFEM library. For more information and source code
-// availability visit https://mfem.org.
-//
-// MFEM is free software; you can redistribute it and/or modify it under the
-// terms of the BSD-3 license. We welcome feedback and contributions, see file
-// CONTRIBUTING.md for details.
-
-#include "../config/config.hpp"
-
-#ifdef MFEM_USE_MPI
-
-#include "fem.hpp"
-
-
-namespace mfem
-{
-
-ParPrmBlockNonlinearForm::ParPrmBlockNonlinearForm(Array<ParFiniteElementSpace *> &pf,
-                                                   Array<ParFiniteElementSpace *> &ppf)
-    :PrmBlockNonlinearForm()
-{
-    pBlockGrad = nullptr;
-    SetParSpaces(pf,ppf);
-}
-
-void ParPrmBlockNonlinearForm::SetParSpaces(Array<ParFiniteElementSpace *> &pf,
-                                            Array<ParFiniteElementSpace *> &pprmf)
-{
-   delete pBlockGrad;
-   pBlockGrad = nullptr;
-
-   for (int s1=0; s1<fes.Size(); ++s1)
-   {
-      for (int s2=0; s2<fes.Size(); ++s2)
-      {
-         delete phBlockGrad(s1,s2);
-      }
-   }
-
-   Array<FiniteElementSpace *> serialSpaces(pf.Size());
-   Array<FiniteElementSpace *> prmserialSpaces(pprmf.Size());
-   for (int s=0; s<pf.Size(); s++)
-   {
-      serialSpaces[s] = (FiniteElementSpace *) pf[s];
-   }
-   for (int s=0; s<pprmf.Size(); s++)
-   {
-      prmserialSpaces[s] = (FiniteElementSpace *) pprmf[s];
-   }
-
-   SetSpaces(serialSpaces,prmserialSpaces);
-
-   phBlockGrad.SetSize(fes.Size(), fes.Size());
-
-   for (int s1=0; s1<fes.Size(); ++s1)
-   {
-      for (int s2=0; s2<fes.Size(); ++s2)
-      {
-         phBlockGrad(s1,s2) = new OperatorHandle(Operator::Hypre_ParCSR);
-      }
-   }
-}
-
-ParFiniteElementSpace * ParPrmBlockNonlinearForm::ParFESpace(int k)
-{
-   return (ParFiniteElementSpace *)fes[k];
-}
-
-const ParFiniteElementSpace *ParPrmBlockNonlinearForm::ParFESpace(int k) const
-{
-   return (const ParFiniteElementSpace *)fes[k];
-}
-
-
-ParFiniteElementSpace * ParPrmBlockNonlinearForm::ParPrmFESpace(int k)
-{
-   return (ParFiniteElementSpace *)prmfes[k];
-}
-
-const ParFiniteElementSpace *ParPrmBlockNonlinearForm::ParPrmFESpace(int k) const
-{
-   return (const ParFiniteElementSpace *)prmfes[k];
-}
-
-// Here, rhs is a true dof vector
-void ParPrmBlockNonlinearForm::SetEssentialBC(const
-                                           Array<Array<int> *>&bdr_attr_is_ess,
-                                           Array<Vector *> &rhs)
-{
-   Array<Vector *> nullarray(fes.Size());
-   nullarray = NULL;
-
-   PrmBlockNonlinearForm::SetEssentialBC(bdr_attr_is_ess, nullarray);
-
-   for (int s = 0; s < fes.Size(); ++s)
-   {
-      if (rhs[s])
-      {
-         rhs[s]->SetSubVector(*ess_tdofs[s], 0.0);
-      }
-   }
-}
-
-void ParPrmBlockNonlinearForm::SetPrmEssentialBC(const
-                                           Array<Array<int> *>&bdr_attr_is_ess,
-                                           Array<Vector *> &rhs)
-{
-    Array<Vector *> nullarray(fes.Size());
-    nullarray = NULL;
-
-    PrmBlockNonlinearForm::SetPrmEssentialBC(bdr_attr_is_ess, nullarray);
-
-    for (int s = 0; s < prmfes.Size(); ++s)
-    {
-       if (rhs[s])
-       {
-          rhs[s]->SetSubVector(*prmess_tdofs[s], 0.0);
-       }
-    }
-}
-
-double ParPrmBlockNonlinearForm::GetEnergy(const Vector &x) const
-{
-   xs_true.Update(x.GetData(), block_trueOffsets);
-   xs.Update(block_offsets);
-
-   for (int s = 0; s < fes.Size(); ++s)
-   {
-      fes[s]->GetProlongationMatrix()->Mult(xs_true.GetBlock(s), xs.GetBlock(s));
-   }
-
-   double enloc = PrmBlockNonlinearForm::GetEnergyBlocked(xs,xdv);
-   double englo = 0.0;
-
-   MPI_Allreduce(&enloc, &englo, 1, MPI_DOUBLE, MPI_SUM,
-                 ParFESpace(0)->GetComm());
-
-   return englo;
-}
-
-void ParPrmBlockNonlinearForm::Mult(const Vector &x, Vector &y) const
-{
-   xs_true.Update(x.GetData(), block_trueOffsets);
-   ys_true.Update(y.GetData(), block_trueOffsets);
-   xs.Update(block_offsets);
-   ys.Update(block_offsets);
-
-   for (int s=0; s<fes.Size(); ++s)
-   {
-      fes[s]->GetProlongationMatrix()->Mult(
-         xs_true.GetBlock(s), xs.GetBlock(s));
-   }
-
-   PrmBlockNonlinearForm::MultBlocked(xs, xdv, ys);
-
-   if (fnfi.Size() > 0)
-   {
-      MFEM_ABORT("TODO: assemble contributions from shared face terms");
-   }
-
-   for (int s=0; s<fes.Size(); ++s)
-   {
-      fes[s]->GetProlongationMatrix()->MultTranspose(
-         ys.GetBlock(s), ys_true.GetBlock(s));
-
-      ys_true.GetBlock(s).SetSubVector(*ess_tdofs[s], 0.0);
-   }
-}
-
-/// Block T-Vector to Block T-Vector
-void ParPrmBlockNonlinearForm::PrmMult(const Vector &x, Vector &y) const
-{
-    xs_true.Update(x.GetData(), prmblock_trueOffsets);
-    ys_true.Update(y.GetData(), prmblock_trueOffsets);
-    prmxs.Update(prmblock_offsets);
-    prmys.Update(prmblock_offsets);
-
-    for (int s=0; s<prmfes.Size(); ++s)
-    {
-       prmfes[s]->GetProlongationMatrix()->Mult(
-          xs_true.GetBlock(s), prmxs.GetBlock(s));
-    }
-
-    PrmBlockNonlinearForm::MultPrmBlocked(xsv,adv,xdv,prmys);
-
-    if (fnfi.Size() > 0)
-    {
-       MFEM_ABORT("TODO: assemble contributions from shared face terms");
-    }
-
-    for (int s=0; s<prmfes.Size(); ++s)
-    {
-       prmfes[s]->GetProlongationMatrix()->MultTranspose(
-          prmys.GetBlock(s), ys_true.GetBlock(s));
-
-       ys_true.GetBlock(s).SetSubVector(*prmess_tdofs[s], 0.0);
-    }
-
-}
-
-/// Return the local gradient matrix for the given true-dof vector x
-const BlockOperator & ParPrmBlockNonlinearForm::GetLocalGradient(
-   const Vector &x) const
-{
-   xs_true.Update(x.GetData(), block_trueOffsets);
-   xs.Update(block_offsets);
-
-   for (int s=0; s<fes.Size(); ++s)
-   {
-      fes[s]->GetProlongationMatrix()->Mult(
-         xs_true.GetBlock(s), xs.GetBlock(s));
-   }
-
-   PrmBlockNonlinearForm::ComputeGradientBlocked(xs,xdv); // (re)assemble Grad with b.c.
-
-   delete BlockGrad;
-   BlockGrad = new BlockOperator(block_offsets);
-
-   for (int i = 0; i < fes.Size(); ++i)
-   {
-      for (int j = 0; j < fes.Size(); ++j)
-      {
-         BlockGrad->SetBlock(i, j, Grads(i, j));
-      }
-   }
-   return *BlockGrad;
-}
-
-// Set the operator type id for the parallel gradient matrix/operator.
-void ParPrmBlockNonlinearForm::SetGradientType(Operator::Type tid)
-{
-   for (int s1=0; s1<fes.Size(); ++s1)
-   {
-      for (int s2=0; s2<fes.Size(); ++s2)
-      {
-         phBlockGrad(s1,s2)->SetType(tid);
-      }
-   }
-}
-
-BlockOperator & ParPrmBlockNonlinearForm::GetGradient(const Vector &x) const
-{
-   if (pBlockGrad == NULL)
-   {
-      pBlockGrad = new BlockOperator(block_trueOffsets);
-   }
-
-   Array<const ParFiniteElementSpace *> pfes(fes.Size());
-
-   for (int s1=0; s1<fes.Size(); ++s1)
-   {
-      pfes[s1] = ParFESpace(s1);
-
-      for (int s2=0; s2<fes.Size(); ++s2)
-      {
-         phBlockGrad(s1,s2)->Clear();
-      }
-   }
-
-   GetLocalGradient(x); // gradients are stored in 'Grads'
-
-   if (fnfi.Size() > 0)
-   {
-      MFEM_ABORT("TODO: assemble contributions from shared face terms");
-   }
-
-   for (int s1=0; s1<fes.Size(); ++s1)
-   {
-      for (int s2=0; s2<fes.Size(); ++s2)
-      {
-         OperatorHandle dA(phBlockGrad(s1,s2)->Type()),
-                        Ph(phBlockGrad(s1,s2)->Type()),
-                        Rh(phBlockGrad(s1,s2)->Type());
-
-         if (s1 == s2)
-         {
-            dA.MakeSquareBlockDiag(pfes[s1]->GetComm(), pfes[s1]->GlobalVSize(),
-                                   pfes[s1]->GetDofOffsets(), Grads(s1,s1));
-            Ph.ConvertFrom(pfes[s1]->Dof_TrueDof_Matrix());
-            phBlockGrad(s1,s1)->MakePtAP(dA, Ph);
-
-            OperatorHandle Ae;
-            Ae.EliminateRowsCols(*phBlockGrad(s1,s1), *ess_tdofs[s1]);
-         }
-         else
-         {
-            dA.MakeRectangularBlockDiag(pfes[s1]->GetComm(),
-                                        pfes[s1]->GlobalVSize(),
-                                        pfes[s2]->GlobalVSize(),
-                                        pfes[s1]->GetDofOffsets(),
-                                        pfes[s2]->GetDofOffsets(),
-                                        Grads(s1,s2));
-            Rh.ConvertFrom(pfes[s1]->Dof_TrueDof_Matrix());
-            Ph.ConvertFrom(pfes[s2]->Dof_TrueDof_Matrix());
-
-            phBlockGrad(s1,s2)->MakeRAP(Rh, dA, Ph);
-
-            phBlockGrad(s1,s2)->EliminateRows(*ess_tdofs[s1]);
-            phBlockGrad(s1,s2)->EliminateCols(*ess_tdofs[s2]);
-         }
-
-         pBlockGrad->SetBlock(s1, s2, phBlockGrad(s1,s2)->Ptr());
-      }
-   }
-
-   return *pBlockGrad;
-}
-
-ParPrmBlockNonlinearForm::~ParPrmBlockNonlinearForm()
-{
-   delete pBlockGrad;
-   for (int s1=0; s1<fes.Size(); ++s1)
-   {
-      for (int s2=0; s2<fes.Size(); ++s2)
-      {
-         delete phBlockGrad(s1,s2);
-      }
-   }
-}
-
-
-void ParPrmBlockNonlinearForm::SetStateFields(const Vector &xv) const
-{
-   xs_true.Update(xv.GetData(), block_trueOffsets);
-   xsv.Update(block_offsets);
-   for (int s=0; s<fes.Size(); ++s)
-   {
-      fes[s]->GetProlongationMatrix()->Mult(
-         xs_true.GetBlock(s), xsv.GetBlock(s));
-   }
-}
-
-
-void ParPrmBlockNonlinearForm::SetAdjointFields(const Vector &av) const
-{
-    xs_true.Update(av.GetData(), block_trueOffsets);
-    adv.Update(block_offsets);
-    for (int s=0; s<fes.Size(); ++s)
-    {
-       fes[s]->GetProlongationMatrix()->Mult(
-          xs_true.GetBlock(s), adv.GetBlock(s));
-    }
-}
-
-void ParPrmBlockNonlinearForm::SetPrmFields(const Vector &dv) const
-{
-    xs_true.Update(dv.GetData(),prmblock_trueOffsets);
-    xdv.Update(prmblock_offsets);
-    for (int s=0; s<prmfes.Size(); ++s)
-    {
-       prmfes[s]->GetProlongationMatrix()->Mult(
-          xs_true.GetBlock(s), xdv.GetBlock(s));
-    }
-}
-
-
-}
-
-#endif

fem/pprmnonlinearform.hpp

@@ -1,105 +0,0 @@
-// Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
-// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
-// LICENSE and NOTICE for details. LLNL-CODE-806117.
-//
-// This file is part of the MFEM library. For more information and source code
-// availability visit https://mfem.org.
-//
-// MFEM is free software; you can redistribute it and/or modify it under the
-// terms of the BSD-3 license. We welcome feedback and contributions, see file
-// CONTRIBUTING.md for details.
-
-#ifndef MFEM_PPRMNONLINEARFORM
-#define MFEM_PPRMNONLINEARFORM
-
-#include "../config/config.hpp"
-
-#ifdef MFEM_USE_MPI
-
-#include "pgridfunc.hpp"
-#include "prmnonlinearform.hpp"
-
-namespace mfem
-{
-
-/** @brief A class representing a general parametric parallel block nonlinear operator
-    defined on the Cartesian product of multiple ParFiniteElementSpace%s. */
-/** The ParPrmBlockNonlinearForm takes as input, and returns as output, vectors on
-    the true dofs. */
-class ParPrmBlockNonlinearForm : public PrmBlockNonlinearForm
-{
-protected:
-   mutable BlockVector xs_true, ys_true;
-   mutable Array2D<OperatorHandle *> phBlockGrad;
-   mutable BlockOperator *pBlockGrad;
-
-public:
-   /// Computes the energy of the system
-   virtual double GetEnergy(const Vector &x) const;
-
-   /// Construct an empty ParPrmBlockNonlinearForm. Initialize with SetParSpaces().
-   ParPrmBlockNonlinearForm() : pBlockGrad(nullptr) { }
-
-   /** @brief Construct a ParPrmBlockNonlinearForm on the given set of
-       parametric and state ParFiniteElementSpace%s. */
-   ParPrmBlockNonlinearForm(Array<ParFiniteElementSpace *> &pf, Array<ParFiniteElementSpace *> &ppf );
-
-   /// Return the @a k-th parallel FE state space of the ParPrmBlockNonlinearForm.
-   ParFiniteElementSpace *ParFESpace(int k);
-   /** @brief Return the @a k-th parallel FE state space of the ParPrmBlockNonlinearForm
-       (const version). */
-   const ParFiniteElementSpace *ParFESpace(int k) const;
-
-   /// Return the @a k-th parallel FE parameters space of the ParPrmBlockNonlinearForm.
-   ParFiniteElementSpace *ParPrmFESpace(int k);
-   /** @brief Return the @a k-th parallel FE parameters space of the ParPrmBlockNonlinearForm
-       (const version). */
-   const ParFiniteElementSpace *ParPrmFESpace(int k) const;
-
-
-   /** @brief After a call to SetParSpaces(), the essential b.c. and the
-       gradient-type (if different from the default) must be set again. */
-   void SetParSpaces(Array<ParFiniteElementSpace *> &pf, Array<ParFiniteElementSpace *> &pprmf);
-
-   // Here, rhs is a true dof vector
-   virtual void SetEssentialBC(const Array<Array<int> *>&bdr_attr_is_ess,
-                                  Array<Vector *> &rhs);
-
-   // Here, rhs is a true dof vector
-   virtual void SetPrmEssentialBC(const Array<Array<int> *>&bdr_attr_is_ess,
-                                  Array<Vector *> &rhs);
-
-
-   /// Block T-Vector to Block T-Vector
-   virtual void Mult(const Vector &x, Vector &y) const;
-
-   /// Block T-Vector to Block T-Vector
-   virtual void PrmMult(const Vector &x, Vector &y) const;
-
-   /// Return the local block gradient matrix for the given true-dof vector x
-   const BlockOperator &GetLocalGradient(const Vector &x) const;
-
-   virtual BlockOperator &GetGradient(const Vector &x) const;
-
-   /** @brief Set the operator type id for the blocks of the parallel gradient
-       matrix/operator. The default type is Operator::Hypre_ParCSR. */
-   void SetGradientType(Operator::Type tid);
-
-   /// Destructor.
-   virtual ~ParPrmBlockNonlinearForm();
-
-   /// Set the state fields
-   virtual void SetStateFields(const Vector &xv) const;
-
-   /// Set the adjoint fields
-   virtual void SetAdjointFields(const Vector &av) const;
-
-   /// Set the parameters/design fields
-   virtual void SetPrmFields(const Vector &dv) const;
-
-};
-
-}
-
-#endif
-#endif

fem/prmnonlinearform.hpp

@@ -1,231 +0,0 @@
-// Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
-// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
-// LICENSE and NOTICE for details. LLNL-CODE-806117.
-//
-// This file is part of the MFEM library. For more information and source code
-// availability visit https://mfem.org.
-//
-// MFEM is free software; you can redistribute it and/or modify it under the
-// terms of the BSD-3 license. We welcome feedback and contributions, see file
-// CONTRIBUTING.md for details.
-
-#ifndef MFEM_PRMNONLINEARFORM
-#define MFEM_PRMNONLINEARFORM
-
-#include "../config/config.hpp"
-#include "nonlininteg.hpp"
-#include "nonlinearform_ext.hpp"
-#include "bilinearform.hpp"
-#include "gridfunc.hpp"
-
-namespace mfem
-{
-
-/** @brief A class representing a general parametric block nonlinear operator
-    defined on the Cartesian product of multiple FiniteElementSpace%s. */
-class PrmBlockNonlinearForm : public Operator
-{
-protected:
-    /// FE spaces on which the form lives.
-    Array<FiniteElementSpace*> fes;
-
-    /// FE spaces for the parametric fields
-    Array<FiniteElementSpace*> prmfes;
-
-    int prmheight;
-    int prmwidth;
-
-
-    /// Set of Domain Integrators to be assembled (added).
-    Array<PrmBlockNonlinearFormIntegrator*> dnfi;
-
-    /// Set of interior face Integrators to be assembled (added).
-    Array<PrmBlockNonlinearFormIntegrator*> fnfi;
-
-    /// Set of Boundary Face Integrators to be assembled (added).
-    Array<PrmBlockNonlinearFormIntegrator*> bfnfi;
-    Array<Array<int>*>           bfnfi_marker;
-
-    /** Auxiliary block-vectors for wrapping input and output vectors or holding
-           GridFunction-like block-vector data (e.g. in parallel). */
-    mutable BlockVector xs, ys;
-    mutable BlockVector prmxs, prmys;
-
-
-    /** Auxiliary block-vectors for holding
-           GridFunction-like block-vector data (e.g. in parallel). */
-    mutable BlockVector xsv;
-
-    /** Auxiliary block-vectors for holding
-           GridFunction-like block-vector data for the parameter fields
-                                                (e.g. in parallel). */
-    mutable BlockVector xdv;
-    /** Auxiliary block-vectors for holding
-           GridFunction-like block-vector data for the adjoint fields
-                                                (e.g. in parallel). */
-    mutable BlockVector adv;
-
-
-    mutable Array2D<SparseMatrix*> Grads, cGrads;
-    mutable BlockOperator *BlockGrad;
-
-    // A list of the offsets
-    Array<int> block_offsets;
-    Array<int> block_trueOffsets;
-    // A list with the offsets for the parametric fields
-    Array<int> prmblock_offsets;
-    Array<int> prmblock_trueOffsets;
-
-    // Array of Arrays of tdofs for each space in 'fes'
-    Array<Array<int> *> ess_tdofs;
-
-    // Array of Arrays of tdofs for each space in 'prmfes'
-    Array<Array<int> *> prmess_tdofs;
-
-    /// Array of pointers to the prolongation matrix of fes, may be NULL
-    Array<const Operator *> P;
-
-    /// Array of pointers to the prolongation matrix of prmfes, may be NULL
-    Array<const Operator *> Pprm;
-
-    /// Array of results of dynamic-casting P to SparseMatrix pointer
-    Array<const SparseMatrix *> cP;
-
-    /// Array of results of dynamic-casting Pprm to SparseMatrix pointer
-    Array<const SparseMatrix *> cPprm;
-
-
-    /// Indicator if the Operator is part of a parallel run
-    bool is_serial = true;
-
-    /// Indicator if the Operator needs prolongation on assembly
-    bool needs_prolongation = false;
-
-    /// Indicator if the Operator needs prolongation on assembly
-    bool prmneeds_prolongation = false;
-
-
-    mutable BlockVector aux1, aux2;
-
-    mutable BlockVector prmaux1, prmaux2;
-
-    const BlockVector &Prolongate(const BlockVector &bx) const;
-
-    const BlockVector &PrmProlongate(const BlockVector &bx) const;
-
-    /// Specialized version of GetEnergy() for BlockVectors
-    //double GetEnergyBlocked(const BlockVector &bx) const;
-    double GetEnergyBlocked(const BlockVector &bx, const BlockVector &dx) const;
-
-
-    /// Specialized version of Mult() for BlockVector%s
-    /// Block L-Vector to Block L-Vector
-    void MultBlocked(const BlockVector &bx, const BlockVector &dx, BlockVector &by) const;
-
-    /// Specialized version of Mult() for BlockVector%s
-    /// Block L-Vector to Block L-Vector
-    /// bx - state vector, ax - adjoint vector, dx - parametric fields
-    /// dy = ax' d(residual(bx))/d(dx)
-    void MultPrmBlocked(const BlockVector &bx, const BlockVector & ax, const BlockVector &dx, BlockVector &dy) const;
-
-
-    /// Specialized version of GetGradient() for BlockVector
-    //void ComputeGradientBlocked(const BlockVector &bx) const;
-    void ComputeGradientBlocked(const BlockVector &bx, const BlockVector &dx) const;
-
-public:
-   /// Construct an empty BlockNonlinearForm. Initialize with SetSpaces().
-   PrmBlockNonlinearForm();
-
-   /// Construct a BlockNonlinearForm on the given set of FiniteElementSpace%s.
-   PrmBlockNonlinearForm(Array<FiniteElementSpace *> &f, Array<FiniteElementSpace *> &pf );
-
-   /// Return the @a k-th FE space of the PrmBlockNonlinearForm.
-   FiniteElementSpace *FESpace(int k) { return fes[k]; }
-
-   /// Return the @a k-th parametric FE space of the PrmBlockNonlinearForm.
-   FiniteElementSpace *PrmFESpace(int k) { return prmfes[k]; }
-
-
-   /// Return the @a k-th FE space of the BlockNonlinearForm (const version).
-   const FiniteElementSpace *FESpace(int k) const { return fes[k]; }
-
-   /// Return the @a k-th parametric FE space of the BlockNonlinearForm (const version).
-   const FiniteElementSpace *PrmFESpace(int k) const { return prmfes[k]; }
-
-   Array<PrmBlockNonlinearFormIntegrator*>& GetDNFI(){ return dnfi;}
-
-
-   /// (Re)initialize the PrmBlockNonlinearForm.
-   /** After a call to SetSpaces(), the essential b.c. must be set again. */
-   void SetSpaces(Array<FiniteElementSpace *> &f, Array<FiniteElementSpace *> &prmf);
-
-   /// Return the regular dof offsets.
-   const Array<int> &GetBlockOffsets() const { return block_offsets; }
-   /// Return the true-dof offsets.
-   const Array<int> &GetBlockTrueOffsets() const { return block_trueOffsets; }
-
-   /// Return the regular dof offsets for the parameters.
-   const Array<int> &PrmGetBlockOffsets() const { return prmblock_offsets; }
-   /// Return the true-dof offsets for the parameters.
-   const Array<int> &PrmGetBlockTrueOffsets() const { return prmblock_trueOffsets; }
-
-   /// Adds new Domain Integrator.
-   void AddDomainIntegrator(PrmBlockNonlinearFormIntegrator *nlfi)
-   { dnfi.Append(nlfi); }
-
-   /// Adds new Interior Face Integrator.
-   void AddInteriorFaceIntegrator(PrmBlockNonlinearFormIntegrator *nlfi)
-   { fnfi.Append(nlfi); }
-
-   /// Adds new Boundary Face Integrator.
-   void AddBdrFaceIntegrator(PrmBlockNonlinearFormIntegrator *nlfi)
-   { bfnfi.Append(nlfi); bfnfi_marker.Append(NULL); }
-
-   /** @brief Adds new Boundary Face Integrator, restricted to specific boundary
-       attributes. */
-   void AddBdrFaceIntegrator(PrmBlockNonlinearFormIntegrator *nlfi,
-                             Array<int> &bdr_marker);
-
-   virtual void SetEssentialBC(const Array<Array<int> *>&bdr_attr_is_ess,
-                               Array<Vector *> &rhs);
-
-   virtual void SetPrmEssentialBC(const Array<Array<int> *>&bdr_attr_is_ess,
-                               Array<Vector *> &rhs);
-
-
-   virtual double GetEnergy(const Vector &x) const;
-
-   /// Method is only called in serial, the parallel version calls MultBlocked
-   /// directly.
-   virtual void Mult(const Vector &x, Vector &y) const;
-
-   /// Method is only called in serial, the parallel version calls MultBlocked
-   /// directly.
-   virtual void PrmMult(const Vector &x, Vector &t) const;
-
-   /// Method is only called in serial, the parallel version calls
-   /// GetGradientBlocked directly.
-   virtual Operator &GetGradient(const Vector &x) const;
-
-   /// Set the state fields
-   virtual void SetStateFields(const Vector &xv) const;
-
-   /// Set the adjoint fields
-   virtual void SetAdjointFields(const Vector &av) const;
-
-   /// Set the parameters/design fields
-   virtual void SetPrmFields(const Vector &dv) const;
-
-
-
-   /// Destructor.
-   virtual ~PrmBlockNonlinearForm();
-
-};
-
-}
-
-#endif
-
-

general/CMakeLists.txt

@@ -27,7 +27,6 @@ list(APPEND SRCS
   stable3d.cpp
   table.cpp
   tic_toc.cpp
-  tinyxml2.cpp
   version.cpp
   )
 
@@ -56,7 +55,6 @@ list(APPEND HDRS
   table.hpp
   tassign.hpp
   tic_toc.hpp
-  tinyxml2.h
   text.hpp
   version.hpp
   )

linalg/CMakeLists.txt

@@ -10,6 +10,7 @@
 # CONTRIBUTING.md for details.
 
 list(APPEND SRCS
+  auxiliary.cpp
   blockmatrix.cpp
   blockoperator.cpp
   blockvector.cpp
@@ -27,6 +28,7 @@ list(APPEND SRCS
   )
 
 list(APPEND HDRS
+  auxiliary.hpp
   blockmatrix.hpp
   blockoperator.hpp
   blockvector.hpp

linalg/auxiliary.hpp

@@ -0,0 +1,285 @@
+// Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef MFEM_AUXILIARY
+#define MFEM_AUXILIARY
+
+#include "../config/config.hpp"
+
+#ifdef MFEM_USE_MPI
+
+#include "../general/tic_toc.hpp"
+#include "solvers.hpp"
+
+namespace mfem
+{
+
+// forward declarations
+class Coefficient;
+class MatrixCoefficient;
+class ParMesh;
+class ParBilinearForm;
+class ParDiscreteLinearOperator;
+
+/** @brief Auxiliary space solvers for MatrixFreeAMS preconditioner
+
+    Given an operator A and a transfer G, this will create a solver
+    that approximates (G^T A G)^{-1}. Used for two different
+    auxiliary spaces in the AMS cycle.
+
+    The produced solver is based on a low-order refined discretization
+    for the high-order H1 problem. */
+class MatrixFreeAuxiliarySpace : public Solver
+{
+public:
+   /** @brief Pi space constructor
+
+       In the AMS framework this auxiliary space has two coefficients.
+
+       @param mesh_lor       Low-order refined auxiliary mesh
+       @param alpha_coeff    coefficient on curl-curl term (1 if null)
+       @param beta_coeff     coefficient on mass term (1 if null)
+       @param beta_mcoeff    matrix coefficient on mass term
+       @param ess_bdr        attributes for essential boundaries
+       @param curlcurl_oper  High-order operator for the system
+       @param pi             Intentity interpolation operator
+       @param useAmgX_       Use AmgX instead of hypre for auxiliary solves
+       @param cg_iterations  number of CG iterations used to invert
+                             auxiliary system, choosing 0 means to use a
+                             single V-cycle
+   */
+   MatrixFreeAuxiliarySpace(
+      ParMesh& mesh_lor, Coefficient* alpha_coeff, Coefficient* beta_coeff,
+      MatrixCoefficient* beta_mcoeff,
+      Array<int>& ess_bdr, Operator& curlcurl_oper, Operator& pi,
+#ifdef MFEM_USE_AMGX
+      bool useAmgX_,
+#endif
+      int cg_iterations = 0);
+
+   // Complex Pi space constructor
+   MatrixFreeAuxiliarySpace(
+      ParMesh& mesh_lor, Coefficient* alpha_coeff, Coefficient* beta_coeff,
+      Coefficient* beta_imag, Coefficient* abs_beta_imag,
+      MatrixCoefficient* beta_mcoeff,
+      Array<int>& ess_bdr, Operator& curlcurl_oper, Operator *oper_complex,
+      Operator& pi,
+#ifdef MFEM_USE_AMGX
+      bool useAmgX_,
+#endif
+      int cg_iterations = 0);
+
+   /** @brief G space constructor
+
+       This has one coefficient in the AMS framework.
+
+       @param mesh_lor       Low-order refined auxiliary mesh
+       @param beta_coeff     coefficient on mass term (1 if null)
+       @param beta_mcoeff    matrix coefficient on mass term
+       @param ess_bdr        attributes for essential boundaries
+       @param curlcurl_oper  High-order operator for the system
+       @param g              Gradient interpolation operator
+       @param useAmgX_       Use AmgX instead of hypre for auxiliary solves
+       @param cg_iterations  number of CG iterations used to invert
+                             auxiliary system, choosing 0 means to
+                             use a single V-cycle
+   */
+   MatrixFreeAuxiliarySpace(
+      ParMesh& mesh_lor, Coefficient* beta_coeff,
+      MatrixCoefficient* beta_mcoeff, Array<int>& ess_bdr,
+      Operator& curlcurl_oper, Operator& g,
+#ifdef MFEM_USE_AMGX
+      bool useAmgX_,
+#endif
+      int cg_iterations = 1);
+
+   // Complex G space constructor
+   MatrixFreeAuxiliarySpace(
+      ParMesh& mesh_lor, Coefficient* beta_coeff, Coefficient* beta_imag,
+      Coefficient* abs_beta_imag,
+      MatrixCoefficient* beta_mcoeff, Array<int>& ess_bdr,
+      Operator& curlcurl_oper, Operator *oper_complex, Operator& g,
+#ifdef MFEM_USE_AMGX
+      bool useAmgX_,
+#endif
+      int cg_iterations = 1);
+
+   ~MatrixFreeAuxiliarySpace();
+
+   void Mult(const Vector& x, Vector& y) const;
+
+   void SetOperator(const Operator& op) {}
+
+private:
+   /** @brief Helper routine for constructors.
+
+       @param system_dimension is passed to HypreBoomerAMG::SetSystemsOptions
+   */
+   void SetupAMG(int system_dimension);
+   void SetupVCycle();
+
+   /// inner_cg_iterations > 99 applies an exact solve here
+   void SetupCG(Operator& curlcurl_oper, Operator& conn,
+                int inner_cg_iterations);
+
+   void SetupGMRES(Operator& curlcurl_oper, Operator& conn);
+
+   void SetupPMHSS();
+
+   MPI_Comm comm;
+   Array<int> ess_tdof_list;
+   HypreParMatrix * aspacematrix;
+   HypreParMatrix * aspacematrix_complex;
+   HypreParMatrix * aspacematrix_imag;
+   Solver * aspacepc;
+   Operator* matfree;
+   CGSolver* cg;
+   GMRESSolver* gmres;
+   GMRESSolver* gmres_PMHSS;
+   Operator* aspacewrapper;
+#ifdef MFEM_USE_AMGX
+   const bool useAmgX;
+#endif
+   mutable int inner_aux_iterations;
+
+   const bool imagBdry;
+
+   Complex_PMHSS *PMHSS = NULL;
+
+   Array<int> offsets;
+   Array<int> offsets_nd;
+   BlockDiagonalPreconditioner *BlockDP;
+
+   BlockOperator *conn_block;
+};
+
+
+/** @brief Perform AMS cycle with generic Operator objects.
+
+    Most users should use MatrixFreeAMS, which wraps this. */
+class GeneralAMS : public Solver
+{
+public:
+   /** @brief Constructor.
+
+       Most of these arguments just need a Mult() operation,
+       but pi and g also require MultTranspose() */
+   GeneralAMS(const Operator& curlcurl_op_,
+              Operator *oper_complex,
+              const Operator& pi_,
+              const Operator& gradient_,
+              const Operator& pispacesolver_,
+              const Operator& gspacesolver_,
+              const Operator& smoother_,
+              const Array<int>& ess_tdof_list_);
+   virtual ~GeneralAMS();
+
+   /// in principle this should set A_ = op;
+   void SetOperator(const Operator &op) {}
+
+   virtual void Mult(const Vector& x, Vector& y) const;
+
+private:
+   const Operator& curlcurl_op;
+   Operator *oper_complex;
+   const Operator& pi;
+   const Operator& gradient;
+   const Operator& pispacesolver;
+   const Operator& gspacesolver;
+   const Operator& smoother;
+   const Array<int> ess_tdof_list;
+
+   void FormResidual(const Vector& rhs, const Vector& x,
+                     Vector& residual) const;
+};
+
+
+/** @brief An auxiliary Maxwell solver for a high-order curl-curl
+    system without high-order assembly.
+
+    The auxiliary space solves are done using a low-order refined approach,
+    but all the interpolation operators, residuals, etc. are done in a
+    matrix-free manner.
+
+    See Barker and Kolev, Matrix-free preconditioning for high-order H(curl)
+    discretizations (https://doi.org/10.1002/nla.2348) */
+class MatrixFreeAMS : public Solver
+{
+public:
+   /** @brief Construct matrix-free AMS preconditioner
+
+       @param aform        BilinearForm for curl-curl problem, generally will
+                           have a CurlCurlIntegrator and possibly a
+                           VectorFEMassIntegrator.
+       @param oper         Operator to precondition.
+       @param nd_fespace   Underlying Nedelec finite element space.
+       @param alpha_coeff  coefficient on curl-curl term in Maxwell problem
+                           (can be null, in which case constant 1 is assumed)
+       @param beta_coeff   (scalar) coefficient on mass term in Maxwell problem
+       @param beta_mcoeff  (matrix) coefficient on mass term
+       @param ess_bdr      boundary *attributes* that are marked essential. In
+                           contrast to other MFEM cases, these are *attributes*
+                           not dofs, because we need to apply these boundary
+                           conditions to different bilinear forms.
+       @param useAmgX      use AmgX (instead of hypre) for LOR problems
+       @param inner_pi_its number of CG iterations on auxiliary pi space,
+                           may need more for difficult coefficients
+       @param inner_g_its  number of CG iterations on auxiliary g space,
+                           may need more for difficult coefficients
+       @param nd_smoother  optional user-provided smoother for Nedelec space,
+                           this object takes ownership and will delete.
+   */
+   MatrixFreeAMS(ParBilinearForm& aform, Operator& oper, Operator *oper_complex,
+                 ParFiniteElementSpace& nd_fespace, Coefficient* alpha_coeff,
+                 Coefficient* beta_coeff, Coefficient* beta_imag,
+                 Coefficient* abs_beta_imag, MatrixCoefficient* beta_mcoeff,
+                 Array<int>& ess_bdr,
+#ifdef MFEM_USE_AMGX
+                 bool useAmgX = false,
+#endif
+                 int inner_pi_its = 0, int inner_g_its = 1,
+                 Solver* nd_smoother = NULL);
+
+   ~MatrixFreeAMS();
+
+   void SetOperator(const Operator &op) {}
+
+   void Mult(const Vector& x, Vector& y) const { general_ams->Mult(x, y); }
+
+private:
+   GeneralAMS * general_ams;
+
+   Solver * smoother;
+   ParDiscreteLinearOperator * pa_grad;
+   OperatorPtr Gradient;
+   ParDiscreteLinearOperator * pa_interp;
+   OperatorPtr Pi;
+
+   Solver * Gspacesolver;
+   Solver * Pispacesolver;
+
+   ParFiniteElementSpace * h1_fespace;
+   ParFiniteElementSpace * h1_fespace_d;
+
+   Array<int> offsets_nd;
+   Array<int> offsets_vector;
+   Array<int> offsets_scalar;
+
+   BlockOperator *Pi_block;
+   BlockOperator *Gradient_block;
+   BlockOperator *smoother_block;
+};
+
+} // namespace mfem
+
+#endif // MFEM_USE_MPI
+
+#endif

linalg/hypre.cpp

@@ -128,15 +128,6 @@ HypreParVector::HypreParVector(ParFiniteElementSpace *pfes)
    own_ParVector = 1;
 }
 
-void HypreParVector::WrapHypreParVector(hypre_ParVector *y, bool owner)
-{
-   if (own_ParVector) { hypre_ParVectorDestroy(x); }
-   Destroy();
-   x = y;
-   _SetDataAndSize_();
-   own_ParVector = owner;
-}
-
 Vector * HypreParVector::GlobalVector() const
 {
    hypre_Vector *hv = hypre_ParVectorToVectorAll(*this);
@@ -935,15 +926,6 @@ void HypreParMatrix::GetOffd(SparseMatrix &offd, HYPRE_Int* &cmap) const
    cmap = A->col_map_offd;
 }
 
-void HypreParMatrix::MergeDiagAndOffd(SparseMatrix &merged)
-{
-   SparseMatrix tmp_wrapper;
-   hypre_CSRMatrix *hypre_merged = hypre_MergeDiagAndOffd(A);
-   MakeWrapper(hypre_merged, tmp_wrapper);
-   merged = tmp_wrapper;
-   hypre_CSRMatrixDestroy(hypre_merged);
-}
-
 void HypreParMatrix::GetBlocks(Array2D<HypreParMatrix*> &blocks,
                                bool interleaved_rows,
                                bool interleaved_cols) const
@@ -984,46 +966,6 @@ HypreParMatrix * HypreParMatrix::Transpose() const
    return new HypreParMatrix(At);
 }
 
-#if MFEM_HYPRE_VERSION >= 21800
-HypreParMatrix *HypreParMatrix::ExtractSubmatrix(const Array<int> &indices,
-                                                 double threshhold) const
-{
-   if (!(A->comm))
-   {
-      hypre_MatvecCommPkgCreate(A);
-   }
-
-   hypre_ParCSRMatrix *submat;
-
-   // Get number of rows stored on this processor
-   int local_num_vars = hypre_CSRMatrixNumRows(hypre_ParCSRMatrixDiag(A));
-
-   // Form hypre CF-splitting array designating submatrix as F-points (-1)
-   Array<int> CF_marker(local_num_vars);
-   CF_marker = 1;
-   for (int j=0; j<indices.Size(); j++)
-   {
-      if (indices[j] > local_num_vars)
-      {
-         MFEM_WARNING("WARNING : " << indices[j] << " > " << local_num_vars);
-      }
-      CF_marker[indices[j]] = -1;
-   }
-
-   // Construct cpts_global array on hypre matrix structure
-   HYPRE_BigInt *cpts_global;
-   hypre_BoomerAMGCoarseParms(MPI_COMM_WORLD, local_num_vars, 1, NULL,
-                              CF_marker, NULL, &cpts_global);
-
-   // Extract submatrix into *submat
-   hypre_ParCSRMatrixExtractSubmatrixFC(A, CF_marker, cpts_global,
-                                        "FF", &submat, threshhold);
-
-   mfem_hypre_TFree(cpts_global);
-   return new HypreParMatrix(submat);
-}
-#endif
-
 HYPRE_Int HypreParMatrix::Mult(HypreParVector &x, HypreParVector &y,
                                double a, double b)
 {
@@ -1627,36 +1569,6 @@ void HypreParMatrix::Destroy()
    }
 }
 
-#if MFEM_HYPRE_VERSION >= 21800
-
-void BlockInverseScale(const HypreParMatrix *A, HypreParMatrix *C,
-                       const Vector *b, HypreParVector *d,
-                       int blocksize, BlockInverseScaleJob job)
-{
-   if (job == BlockInverseScaleJob::MATRIX_ONLY ||
-       job == BlockInverseScaleJob::MATRIX_AND_RHS)
-   {
-      hypre_ParCSRMatrix *C_hypre;
-      hypre_ParcsrBdiagInvScal(*A, blocksize, &C_hypre);
-      hypre_ParCSRMatrixDropSmallEntries(C_hypre, 1e-15, 1);
-      C->WrapHypreParCSRMatrix(C_hypre);
-   }
-
-   if (job == BlockInverseScaleJob::RHS_ONLY ||
-       job == BlockInverseScaleJob::MATRIX_AND_RHS)
-   {
-      HypreParVector b_Hypre(A->GetComm(),
-                             A->GetGlobalNumRows(),
-                             b->GetData(), A->GetRowStarts());
-      hypre_ParVector *d_hypre;
-      hypre_ParvecBdiagInvScal(b_Hypre, blocksize, &d_hypre, *A);
-
-      d->WrapHypreParVector(d_hypre, true);
-   }
-}
-
-#endif
-
 #if MFEM_HYPRE_VERSION < 21400
 
 HypreParMatrix *Add(double alpha, const HypreParMatrix &A,
@@ -2905,11 +2817,6 @@ void HypreGMRES::SetTol(double tol)
    HYPRE_GMRESSetTol(gmres_solver, tol);
 }
 
-void HypreGMRES::SetAbsTol(double tol)
-{
-   HYPRE_GMRESSetAbsoluteTol(gmres_solver, tol);
-}
-
 void HypreGMRES::SetMaxIter(int max_iter)
 {
    HYPRE_GMRESSetMaxIter(gmres_solver, max_iter);
@@ -3713,116 +3620,6 @@ void HypreBoomerAMG::SetElasticityOptions(ParFiniteElementSpace *fespace)
    error_mode = IGNORE_HYPRE_ERRORS;
 }
 
-#if MFEM_HYPRE_VERSION >= 21800
-
-void HypreBoomerAMG::SetAdvectiveOptions(int distanceR,
-                                         const std::string &prerelax,
-                                         const std::string &postrelax)
-{
-   // Hypre parameters
-   int Sabs = 0;
-   int interp_type = 100;
-   int relax_type = 10;
-   int coarsen_type = 6;
-   double strength_tolC = 0.1;
-   double strength_tolR = 0.01;
-   double filter_tolR = 0.0;
-   double filterA_tol = 0.0;
-
-   // Set relaxation on specified grid points
-   int ns_down, ns_up, ns_coarse;
-   if (distanceR > 0)
-   {
-      ns_down = prerelax.length();
-      ns_up = postrelax.length();
-      ns_coarse = 1;
-
-      // Array to store relaxation scheme and pass to Hypre
-      HYPRE_Int **grid_relax_points = mfem_hypre_TAlloc(HYPRE_Int*, 4);
-      grid_relax_points[0] = NULL;
-      grid_relax_points[1] = mfem_hypre_TAlloc(HYPRE_Int, ns_down);
-      grid_relax_points[2] = mfem_hypre_TAlloc(HYPRE_Int, ns_up);
-      grid_relax_points[3] = mfem_hypre_TAlloc(HYPRE_Int, 1);
-      grid_relax_points[3][0] = 0;
-
-      // set down relax scheme
-      for (int i = 0; i<ns_down; i++)
-      {
-         if (prerelax[i] == 'F')
-         {
-            grid_relax_points[1][i] = -1;
-         }
-         else if (prerelax[i] == 'C')
-         {
-            grid_relax_points[1][i] = 1;
-         }
-         else if (prerelax[i] == 'A')
-         {
-            grid_relax_points[1][i] = 0;
-         }
-      }
-
-      // set up relax scheme
-      for (int i = 0; i<ns_up; i++)
-      {
-         if (postrelax[i] == 'F')
-         {
-            grid_relax_points[2][i] = -1;
-         }
-         else if (postrelax[i] == 'C')
-         {
-            grid_relax_points[2][i] = 1;
-         }
-         else if (postrelax[i] == 'A')
-         {
-            grid_relax_points[2][i] = 0;
-         }
-      }
-
-      HYPRE_BoomerAMGSetRestriction(amg_precond, distanceR);
-
-      HYPRE_BoomerAMGSetGridRelaxPoints(amg_precond, grid_relax_points);
-
-      HYPRE_BoomerAMGSetInterpType(amg_precond, interp_type);
-   }
-
-   if (Sabs)
-   {
-      HYPRE_BoomerAMGSetSabs(amg_precond, Sabs);
-   }
-
-   HYPRE_BoomerAMGSetCoarsenType(amg_precond, coarsen_type);
-
-   // does not support aggressive coarsening
-   HYPRE_BoomerAMGSetAggNumLevels(amg_precond, 0);
-
-   HYPRE_BoomerAMGSetStrongThreshold(amg_precond, strength_tolC);
-
-   if (distanceR > 0)
-   {
-      HYPRE_BoomerAMGSetStrongThresholdR(amg_precond, strength_tolR);
-      HYPRE_BoomerAMGSetFilterThresholdR(amg_precond, filter_tolR);
-   }
-
-   if (relax_type > -1)
-   {
-      HYPRE_BoomerAMGSetRelaxType(amg_precond, relax_type);
-   }
-
-   if (distanceR > 0)
-   {
-      HYPRE_BoomerAMGSetCycleNumSweeps(amg_precond, ns_coarse, 3);
-      HYPRE_BoomerAMGSetCycleNumSweeps(amg_precond, ns_down,   1);
-      HYPRE_BoomerAMGSetCycleNumSweeps(amg_precond, ns_up,     2);
-
-      HYPRE_BoomerAMGSetADropTol(amg_precond, filterA_tol);
-      // type = -1: drop based on row inf-norm
-      HYPRE_BoomerAMGSetADropType(amg_precond, -1);
-   }
-}
-
-#endif
-
 HypreBoomerAMG::~HypreBoomerAMG()
 {
    for (int i = 0; i < rbms.Size(); i++)

linalg/hypre.hpp

@@ -81,14 +81,6 @@ private:
    inline void _SetDataAndSize_();
 
 public:
-
-   /// Default constructor, no underlying @a hypre_ParVector is created.
-   HypreParVector()
-   {
-      own_ParVector = false;
-      x = NULL;
-   }
-
    /** @brief Creates vector with given global size and parallel partitioning of
        the rows/columns given by @a col. */
    /** @anchor hypre_partitioning_descr
@@ -121,9 +113,6 @@ public:
    /// MPI communicator
    MPI_Comm GetComm() { return x->comm; }
 
-   /// Converts hypre's format to HypreParVector
-   void WrapHypreParVector(hypre_ParVector *y, bool owner=true);
-
    /// Returns the parallel row/column partitioning
    /** See @ref hypre_partitioning_descr "here" for a description of the
        partitioning array. */
@@ -239,24 +228,15 @@ public:
    /// An empty matrix to be used as a reference to an existing matrix
    HypreParMatrix();
 
-   /// Converts hypre's format to HypreParMatrix
-   /** If @a owner is false, ownership of @a a is not transferred */
-   void WrapHypreParCSRMatrix(hypre_ParCSRMatrix *a, bool owner = true)
-   {
-      Destroy();
-      Init();
-      A = a;
-      ParCSROwner = owner;
-      height = GetNumRows();
-      width = GetNumCols();
-   }
-
    /// Converts hypre's format to HypreParMatrix
    /** If @a owner is false, ownership of @a a is not transferred */
    explicit HypreParMatrix(hypre_ParCSRMatrix *a, bool owner = true)
    {
       Init();
-      WrapHypreParCSRMatrix(a, owner);
+      A = a;
+      if (!owner) { ParCSROwner = 0; }
+      height = GetNumRows();
+      width = GetNumCols();
    }
 
    /// Creates block-diagonal square parallel matrix.
@@ -407,13 +387,6 @@ public:
    void GetDiag(SparseMatrix &diag) const;
    /// Get the local off-diagonal block. NOTE: 'offd' will not own any data.
    void GetOffd(SparseMatrix &offd, HYPRE_Int* &cmap) const;
-   /** @brief Get a single SparseMatrix containing all rows from this processor,
-       merged from the diagonal and off-diagonal blocks stored by the
-       HypreParMatrix. */
-   /** @note The number of columns in the SparseMatrix will be the global number
-       of columns in the parallel matrix, so using this method may result in an
-       integer overflow in the column indices. */
-   void MergeDiagAndOffd(SparseMatrix &merged);
 
    /** Split the matrix into M x N equally sized blocks of parallel matrices.
        The size of 'blocks' must already be set to M x N. */
@@ -424,13 +397,6 @@ public:
    /// Returns the transpose of *this
    HypreParMatrix * Transpose() const;
 
-   /** Returns principle submatrix given by array of indices of connections
-       with relative size > @a threshold in *this. */
-#if MFEM_HYPRE_VERSION >= 21800
-   HypreParMatrix *ExtractSubmatrix(const Array<int> &indices,
-                                    double threshhold=0.0) const;
-#endif
-
    /// Returns the number of rows in the diagonal block of the ParCSRMatrix
    int GetNumRows() const
    {
@@ -583,23 +549,6 @@ public:
    Type GetType() const { return Hypre_ParCSR; }
 };
 
-#if MFEM_HYPRE_VERSION >= 21800
-
-enum class BlockInverseScaleJob
-{
-   MATRIX_ONLY,
-   RHS_ONLY,
-   MATRIX_AND_RHS
-};
-
-/** Constructs and applies block diagonal inverse of HypreParMatrix.
-    The enum @a job specifies whether the matrix or the RHS should be
-    scaled (or both). */
-void BlockInverseScale(const HypreParMatrix *A, HypreParMatrix *C,
-                       const Vector *b, HypreParVector *d,
-                       int blocksize, BlockInverseScaleJob job);
-#endif
-
 /** @brief Return a new matrix `C = alpha*A + beta*B`, assuming that both `A`
     and `B` use the same row and column partitions and the same `col_map_offd`
     arrays. */
@@ -700,12 +649,11 @@ public:
        4    = truncated l1-scaled block Gauss-Seidel/SSOR
        5    = lumped Jacobi
        6    = Gauss-Seidel
-       10   = On-processor forward solve for matrix w/ triangular structure
        16   = Chebyshev
        1001 = Taubin polynomial smoother
        1002 = FIR polynomial smoother. */
    enum Type { Jacobi = 0, l1Jacobi = 1, l1GS = 2, l1GStr = 4, lumpedJacobi = 5,
-               GS = 6, OPFS = 10, Chebyshev = 16, Taubin = 1001, FIR = 1002
+               GS = 6, Chebyshev = 16, Taubin = 1001, FIR = 1002
              };
 
    HypreSmoother();
@@ -819,28 +767,6 @@ public:
    virtual ~HypreSolver();
 };
 
-
-#if MFEM_HYPRE_VERSION >= 21800
-/** Preconditioner for HypreParMatrices that are triangular in some ordering.
-   Finds correct ordering and performs forward substitution on processor
-   as approximate inverse. Exact on one processor. */
-class HypreTriSolve : public HypreSolver
-{
-public:
-   HypreTriSolve() : HypreSolver() { }
-   explicit HypreTriSolve(HypreParMatrix &A) : HypreSolver(&A) { }
-   virtual operator HYPRE_Solver() const { return NULL; }
-
-   virtual HYPRE_PtrToParSolverFcn SetupFcn() const
-   { return (HYPRE_PtrToParSolverFcn) HYPRE_ParCSROnProcTriSetup; }
-   virtual HYPRE_PtrToParSolverFcn SolveFcn() const
-   { return (HYPRE_PtrToParSolverFcn) HYPRE_ParCSROnProcTriSolve; }
-
-   HypreParMatrix* GetData() { return A; }
-   virtual ~HypreTriSolve() { }
-};
-#endif
-
 /// PCG solver in hypre
 class HyprePCG : public HypreSolver
 {
@@ -918,7 +844,6 @@ public:
    virtual void SetOperator(const Operator &op);
 
    void SetTol(double tol);
-   void SetAbsTol(double tol);
    void SetMaxIter(int max_iter);
    void SetKDim(int dim);
    void SetLogging(int logging);
@@ -1203,94 +1128,9 @@ public:
        construct A. */
    void SetElasticityOptions(ParFiniteElementSpace *fespace);
 
-#if MFEM_HYPRE_VERSION >= 21800
-   /** Hypre parameters to use AIR AMG solve for advection-dominated problems.
-       See "Nonsymmetric Algebraic Multigrid Based on Local Approximate Ideal
-       Restriction (AIR)," Manteuffel, Ruge, Southworth, SISC (2018),
-       DOI:/10.1137/17M1144350. Options: "distanceR" -> distance of neighbor
-       DOFs to buld restriction operator; options include 1, 2, and 15 (1.5).
-       Strings "prerelax" and "postrelax" indicate points to relax on:
-       F = F-points, C = C-points, A = all points. E.g., FFC -> relax on
-       F-points, relax again on F-points, then relax on C-points. */
-   void SetAdvectiveOptions(int distance=15,  const std::string &prerelax="",
-                            const std::string &postrelax="FFC");
-
-   /// Expert option - consult hypre documentation/team
-   void SetStrongThresholdR(double strengthR)
-   { HYPRE_BoomerAMGSetStrongThresholdR(amg_precond, strengthR); }
-
-   /// Expert option - consult hypre documentation/team
-   void SetFilterThresholdR(double filterR)
-   { HYPRE_BoomerAMGSetFilterThresholdR(amg_precond, filterR); }
-
-   /// Expert option - consult hypre documentation/team
-   void SetRestriction(int restrict_type)
-   { HYPRE_BoomerAMGSetRestriction(amg_precond, restrict_type); }
-
-   /// Expert option - consult hypre documentation/team
-   void SetIsTriangular()
-   { HYPRE_BoomerAMGSetIsTriangular(amg_precond, 1); }
-
-   /// Expert option - consult hypre documentation/team
-   void SetGMRESSwitchR(int gmres_switch)
-   { HYPRE_BoomerAMGSetGMRESSwitchR(amg_precond, gmres_switch); }
-
-   /// Expert option - consult hypre documentation/team
-   void SetCycleNumSweeps(int prerelax, int postrelax)
-   {
-      HYPRE_BoomerAMGSetCycleNumSweeps(amg_precond, prerelax,  1);
-      HYPRE_BoomerAMGSetCycleNumSweeps(amg_precond, postrelax, 2);
-   }
-#endif
-
    void SetPrintLevel(int print_level)
    { HYPRE_BoomerAMGSetPrintLevel(amg_precond, print_level); }
 
-   void SetMaxIter(int max_iter)
-   { HYPRE_BoomerAMGSetMaxIter(amg_precond, max_iter); }
-
-   /// Expert option - consult hypre documentation/team
-   void SetMaxLevels(int max_levels)
-   { HYPRE_BoomerAMGSetMaxLevels(amg_precond, max_levels); }
-
-   /// Expert option - consult hypre documentation/team
-   void SetTol(double tol)
-   { HYPRE_BoomerAMGSetTol(amg_precond, tol); }
-
-   /// Expert option - consult hypre documentation/team
-   void SetStrengthThresh(double strength)
-   { HYPRE_BoomerAMGSetStrongThreshold(amg_precond, strength); }
-
-   /// Expert option - consult hypre documentation/team
-   void SetInterpolation(int interp_type)
-   { HYPRE_BoomerAMGSetInterpType(amg_precond, interp_type); }
-
-   /// Expert option - consult hypre documentation/team
-   void SetCoarsening(int coarsen_type)
-   { HYPRE_BoomerAMGSetCoarsenType(amg_precond, coarsen_type); }
-
-   /// Expert option - consult hypre documentation/team
-   void SetRelaxType(int relax_type)
-   { HYPRE_BoomerAMGSetRelaxType(amg_precond, relax_type); }
-
-   /// Expert option - consult hypre documentation/team
-   void SetCycleType(int cycle_type)
-   { HYPRE_BoomerAMGSetCycleType(amg_precond, cycle_type); }
-
-   void GetNumIterations(int &num_it)
-   { HYPRE_BoomerAMGGetNumIterations(amg_precond, &num_it); }
-
-   /// Expert option - consult hypre documentation/team
-   void SetNodal(int blocksize)
-   {
-      HYPRE_BoomerAMGSetNumFunctions(amg_precond, blocksize);
-      HYPRE_BoomerAMGSetNodal(amg_precond, 1);
-   }
-
-   /// Expert option - consult hypre documentation/team
-   void SetAggressiveCoarsening(int num_levels)
-   { HYPRE_BoomerAMGSetAggNumLevels(amg_precond, num_levels); }
-
    /// The typecast to HYPRE_Solver returns the internal amg_precond
    virtual operator HYPRE_Solver() const { return amg_precond; }
 
@@ -1299,8 +1139,6 @@ public:
    virtual HYPRE_PtrToParSolverFcn SolveFcn() const
    { return (HYPRE_PtrToParSolverFcn) HYPRE_BoomerAMGSolve; }
 
-   using HypreSolver::Mult;
-
    virtual ~HypreBoomerAMG();
 };
 

linalg/linalg.hpp

@@ -29,6 +29,7 @@
 #include "solvers.hpp"
 #include "handle.hpp"
 #include "invariants.hpp"
+#include "auxiliary.hpp"
 
 #ifdef MFEM_USE_AMGX
 #include "amgxsolver.hpp"

linalg/operator.hpp

@@ -696,6 +696,42 @@ public:
    { A_.Mult(x, y); y *= a_; }
 };
 
+/// General sum operator: x -> A(x)+B(x)
+class SumOperator : public Operator
+{
+   const Operator *A, *B;
+   bool ownA, ownB;
+   mutable Vector z, w;
+   double cA, cB;
+
+public:
+   SumOperator(const Operator *A_, const Operator *B_,
+               bool ownA_, bool ownB_, double cA_, double cB_)
+      : Operator(A_->Height(), B_->Width()),
+        A(A_), B(B_), ownA(ownA_), ownB(ownB_), z(A_->Height()), w(A_->Width()),
+        cA(cA_), cB(cB_)
+   {
+      MFEM_VERIFY(A->Width() == B->Width() && A->Height() == B->Height(),
+                  "incompatible Operators: A->Width() = " << A->Width()
+                  << ", B->Height() = " << B->Height());
+
+      z.UseDevice(true);
+      w.UseDevice(true);
+   }
+
+   ~SumOperator()
+   {
+      if (ownA) { delete A; }
+      if (ownB) { delete B; }
+   }
+
+   virtual void Mult(const Vector &x, Vector &y) const
+   { B->Mult(x, z); A->Mult(x, y); y *= cA; z *= cB; y += z;}
+
+   virtual void MultTranspose(const Vector &x, Vector &y) const
+   { B->MultTranspose(x, w); A->MultTranspose(x, y); y *= cA; w *= cB; y += w;}
+
+};
 
 /** @brief The transpose of a given operator. Switches the roles of the methods
     Mult() and MultTranspose(). */

linalg/solvers.hpp

@@ -852,6 +852,169 @@ public:
 
 #endif // MFEM_USE_SUITESPARSE
 
+class Complex_PMHSS : public Solver
+{
+public:
+   Complex_PMHSS(Operator *Re, Operator *Im, Solver *prec_Re, Solver *prec_Im,
+                 double a_)
+      : Solver(2*Re->Height()), a(a_), A(Re, Im, false, false),
+        A_Re(Re, NULL, false, false),
+        A_Im(Im, NULL, false, false), u(2*Re->Height()), rhs(2*Re->Height()),
+        n(Re->Height())
+   {
+      MFEM_VERIFY(Re->Height() == Im->Height() && Re->Height() == Re->Width() &&
+                  Im->Height() == Im->Width(), "");
+      MFEM_VERIFY(this->Height() == A.Height(), "");
+
+      // Create CG solver for real operator aV + A_Re in complex space.
+
+      V = useIdentityV ? (Operator*) new IdentityOperator(this->Height()) :
+          (Operator*) &A_Re;
+
+      // In the case V = A_Re, it is faster to use a scaled operator than a SumOperator
+      Operator *sumOpRe = useIdentityV ? (Operator*) new SumOperator(V, &A_Re, false,
+                                                                     false, a, 1.0)
+                          : (Operator*) new ScaledOperator(&A_Re, a + 1.0);
+
+      SumOperator *sumOpIm = new SumOperator(V, &A_Im, false, false, a, 1.0);
+
+      CGSolver *cg = new CGSolver(MPI_COMM_WORLD);
+      cg->SetRelTol(1e-6);
+      cg->SetMaxIter(1000);
+      cg->SetPrintLevel(0);
+      cg->SetOperator(*sumOpRe);
+      cg->SetPreconditioner(*prec_Re);
+      cg->iterative_mode = false;
+
+      SRe = cg;
+
+      CGSolver *cgi = new CGSolver(MPI_COMM_WORLD);
+      cgi->SetRelTol(1e-6);
+      cgi->SetMaxIter(1000);
+      cgi->SetPrintLevel(0);
+      cgi->SetOperator(*sumOpIm);
+      if (prec_Im && useIdentityV) { cgi->SetPreconditioner(*prec_Im); }
+      if (!useIdentityV) { cgi->SetPreconditioner(*prec_Re); }
+      cgi->iterative_mode = false;
+
+      /*
+      // For negative definite imaginary part, but then PMHSS does not work?
+      MINRESSolver *cgi = new MINRESSolver(MPI_COMM_WORLD);
+      cgi->SetRelTol(1e-12);
+      cgi->SetMaxIter(1000);
+      cgi->SetPrintLevel(0);
+      cgi->SetOperator(*sumOpIm);
+      if (prec_Im) cgi->SetPreconditioner(*prec_Im);
+      */
+
+      SIm = cgi;
+   }
+
+   void SetOperator(const Operator &op)
+   {
+      MFEM_VERIFY(false, "Don't call SetOperator");
+   }
+
+   void ComputeResidual(const Vector &b, const Vector &sol, Vector &res) const
+   {
+      A.Mult(sol, res);
+      res -= b;
+   }
+
+   void Mult(const Vector &x, Vector &y) const
+   {
+      if (!(x.Size() == Height() && y.Size() == Height()))
+      {
+         std::cout << "bug";
+      }
+
+      MFEM_VERIFY(x.Size() == Height() && y.Size() == Height(), "");
+
+      const double initNorm = x.Norml2();
+      mfem::out << "MHSS RHS norm " << initNorm << '\n';
+
+      // With V = I, use modified HSS (MHSS) from Bai, Benzi, Chen 2010.
+      y = 0.0;
+
+      for (int it=0; it<maxiter; ++it)
+      {
+         // Solve (aI + Re) u = (aI - i Im) y + x
+
+         if (it == 0)
+         {
+            // Optimize the first iteration, when the initial guess is y=0.
+            SRe->Mult(x, u);
+         }
+         else
+         {
+            A_Im.Mult(y, u);  // u = Im y
+            // Set rhs = -i Im y = -i u
+            for (int j=0; j<n; ++j)
+            {
+               rhs[j] = u[n+j];
+               rhs[n+j] = -u[j];
+            }
+
+            rhs += x;
+
+            V->Mult(y, u);
+            rhs.Add(a, u);
+
+            SRe->Mult(rhs, u);
+         }
+
+         // Solve (aI + Im) y = (aI + i Re) u - i x
+
+         A_Re.Mult(u, y);  // y = Re u
+         // Set rhs = i (Re u - x) = i (y - x)
+         for (int j=0; j<n; ++j)
+         {
+            rhs[j] = -(y[n+j] - x[n+j]);
+            rhs[n+j] = y[j] - x[j];
+         }
+
+         if (useIdentityV)
+         {
+            //V->Mult(u, y);
+            //rhs.Add(a, y);
+            rhs.Add(a, u);
+         }
+         else
+         {
+            // Using V = A_Re
+            rhs.Add(a, y);
+         }
+
+         SIm->Mult(rhs, y);
+
+         ComputeResidual(x, y, rhs);
+         const double resNorm = rhs.Norml2();
+         mfem::out << "MHSS iter " << it << " residual norm " << resNorm << '\n';
+
+         if (resNorm / initNorm < tol)
+         {
+            mfem::out << "MHSS converged\n";
+            break;
+         }
+      }
+   }
+
+private:
+   const double a;
+   const int maxiter = 1;
+   ComplexOperator A, A_Re, A_Im;
+   mutable Vector u, rhs;
+   const int n;
+
+   const double tol = 1.0e-8;
+
+   const bool useIdentityV = false;
+   Operator *V = NULL;
+
+   Solver *SRe = NULL;
+   Solver *SIm = NULL;
+};
+
 }
 
 #endif // MFEM_SOLVERS

makefile

@@ -684,10 +684,9 @@ status info:
 ASTYLE_BIN = astyle
 ASTYLE = $(ASTYLE_BIN) --options=$(SRC)config/mfem.astylerc
 ASTYLE_VER = "Artistic Style Version 2.05.1"
-FORMAT_FILES := $(foreach dir,$(DIRS) $(EM_DIRS) config,$(dir)/*.?pp)
-FORMAT_FILES += tests/unit/*.cpp
-FORMAT_FILES += $(foreach dir,general linalg mesh fem,tests/unit/$(dir)/*.?pp)
-FORMAT_FILES := $(filter-out general/tinyxml2.cpp,$(wildcard $(FORMAT_FILES)))
+FORMAT_FILES = $(foreach dir,$(DIRS) $(EM_DIRS) config,"$(dir)/*.?pp")
+FORMAT_FILES += "tests/unit/*.cpp"
+FORMAT_FILES += $(foreach dir,general linalg mesh fem,"tests/unit/$(dir)/*.?pp")
 
 COUT_CERR_FILES = $(foreach dir,$(DIRS),$(dir)/*.[ch]pp)
 COUT_CERR_EXCLUDE = '^general/error\.cpp' '^general/globals\.[ch]pp'

mesh/mesh.cpp

@@ -3448,7 +3448,6 @@ void Mesh::Loader(std::istream &input, int generate_edges,
 
    Clear();
 
-   istream::pos_type beginning_pos = input.tellg();
    string mesh_type;
    input >> ws;
    getline(input, mesh_type);
@@ -3490,20 +3489,11 @@ void Mesh::Loader(std::istream &input, int generate_edges,
    {
       ReadTrueGridMesh(input);
    }
-   else if (mesh_type.rfind("# vtk DataFile Version") == 0)
+   else if (mesh_type == "# vtk DataFile Version 3.0" ||
+            mesh_type == "# vtk DataFile Version 2.0") // VTK
    {
-      int major_vtk_version = mesh_type[mesh_type.length()-3] - '0';
-      // int minor_vtk_version = mesh_type[mesh_type.length()-1] - '0';
-      MFEM_VERIFY(major_vtk_version >= 2 && major_vtk_version <= 4,
-                  "Unsupported VTK format");
       ReadVTKMesh(input, curved, read_gf, finalize_topo);
    }
-   else if (mesh_type.rfind("<VTKFile ") == 0)
-   {
-      // Go back to beginning of stream
-      input.seekg(beginning_pos);
-      ReadXML_VTKMesh(input, curved, read_gf, finalize_topo);
-   }
    else if (mesh_type == "MFEM NURBS mesh v1.0")
    {
       ReadNURBSMesh(input, curved, read_gf);
@@ -8948,7 +8938,7 @@ void Mesh::PrintVTK(std::ostream &out)
          const int nv = elements[i]->GetNVertices();
          out << nv;
          Geometry::Type geom = elements[i]->GetGeometryType();
-         const int *perm = VTKGeometry::VertexPermutation[geom];
+         const int *perm = (geom == Geometry::PRISM) ? vtk_prism_perm : NULL;
          for (int j = 0; j < nv; j++)
          {
             out << ' ' << v[perm ? perm[j] : j];
@@ -9033,9 +9023,35 @@ void Mesh::PrintVTK(std::ostream &out)
    for (int i = 0; i < NumOfElements; i++)
    {
       int vtk_cell_type = 5;
-      Geometry::Type geom = GetElement(i)->GetGeometryType();
-      if (order == 1) { vtk_cell_type = VTKGeometry::Map[geom]; }
-      else if (order == 2) { vtk_cell_type = VTKGeometry::QuadraticMap[geom]; }
+      Geometry::Type geom_type = GetElement(i)->GetGeometryType();
+      if (order == 1)
+      {
+         switch (geom_type)
+         {
+            case Geometry::POINT:        vtk_cell_type = 1;   break;
+            case Geometry::SEGMENT:      vtk_cell_type = 3;   break;
+            case Geometry::TRIANGLE:     vtk_cell_type = 5;   break;
+            case Geometry::SQUARE:       vtk_cell_type = 9;   break;
+            case Geometry::TETRAHEDRON:  vtk_cell_type = 10;  break;
+            case Geometry::CUBE:         vtk_cell_type = 12;  break;
+            case Geometry::PRISM:        vtk_cell_type = 13;  break;
+            default: break;
+         }
+      }
+      else if (order == 2)
+      {
+         switch (geom_type)
+         {
+            case Geometry::SEGMENT:      vtk_cell_type = 21;  break;
+            case Geometry::TRIANGLE:     vtk_cell_type = 22;  break;
+            case Geometry::SQUARE:       vtk_cell_type = 28;  break;
+            case Geometry::TETRAHEDRON:  vtk_cell_type = 24;  break;
+            case Geometry::CUBE:         vtk_cell_type = 29;  break;
+            case Geometry::PRISM:        vtk_cell_type = 32;  break;
+            default: break;
+         }
+      }
+
       out << vtk_cell_type << '\n';
    }
 
@@ -9253,7 +9269,7 @@ void Mesh::PrintVTU(std::ostream &out, int ref, VTKFormat format,
          {
             coff = coff+nv;
             offset.push_back(coff);
-            const int *p = VTKGeometry::VertexPermutation[geom];
+            const int *p = (geom == Geometry::PRISM) ? vtk_prism_perm : NULL;
             for (int k = 0; k < nv; k++, j++)
             {
                WriteBinaryOrASCII(out, buf, np + RG[p ? p[j] : j], " ", format);
@@ -9284,14 +9300,39 @@ void Mesh::PrintVTU(std::ostream &out, int ref, VTKFormat format,
    out << "<DataArray type=\"UInt8\" Name=\"types\" format=\""
        << fmt_str << "\">" << std::endl;
    // cell types
-   const int *vtk_geom_map =
-      high_order_output ? VTKGeometry::HighOrderMap : VTKGeometry::Map;
    for (int i = 0; i < ne; i++)
    {
       Geometry::Type geom = get_geom(i);
       uint8_t vtk_cell_type = 5;
 
-      vtk_cell_type = vtk_geom_map[geom];
+      // VTK element types defined at: https://git.io/JvZLm
+      switch (geom)
+      {
+         case Geometry::POINT:
+            vtk_cell_type = 1;
+            break;
+         case Geometry::SEGMENT:
+            vtk_cell_type = high_order_output ? 68 : 3;
+            break;
+         case Geometry::TRIANGLE:
+            vtk_cell_type = high_order_output ? 69 : 5;
+            break;
+         case Geometry::SQUARE:
+            vtk_cell_type = high_order_output ? 70 : 9;
+            break;
+         case Geometry::TETRAHEDRON:
+            vtk_cell_type = high_order_output ? 71 : 10;
+            break;
+         case Geometry::CUBE:
+            vtk_cell_type = high_order_output ? 72 : 12;
+            break;
+         case Geometry::PRISM:
+            vtk_cell_type = high_order_output ? 73 : 13;
+            break;
+         default:
+            MFEM_ABORT("Unrecognized VTK element type \"" << geom << "\"");
+            break;
+      }
 
       if (high_order_output)
       {
@@ -9440,7 +9481,21 @@ void Mesh::PrintVTK(std::ostream &out, int ref, int field_data)
       int nv = Geometries.GetVertices(geom)->GetNPoints();
       RefG = GlobGeometryRefiner.Refine(geom, ref, 1);
       Array<int> &RG = RefG->RefGeoms;
-      int vtk_cell_type = VTKGeometry::Map[geom];
+      int vtk_cell_type = 5;
+
+      switch (geom)
+      {
+         case Geometry::POINT:        vtk_cell_type = 1;   break;
+         case Geometry::SEGMENT:      vtk_cell_type = 3;   break;
+         case Geometry::TRIANGLE:     vtk_cell_type = 5;   break;
+         case Geometry::SQUARE:       vtk_cell_type = 9;   break;
+         case Geometry::TETRAHEDRON:  vtk_cell_type = 10;  break;
+         case Geometry::CUBE:         vtk_cell_type = 12;  break;
+         case Geometry::PRISM:        vtk_cell_type = 13;  break;
+         default:
+            MFEM_ABORT("Unrecognized VTK element type \"" << geom << "\"");
+            break;
+      }
 
       for (int j = 0; j < RG.Size(); j += nv)
       {

mesh/mesh.hpp

@@ -239,15 +239,8 @@ protected:
    void ReadNetgen2DMesh(std::istream &input, int &curved);
    void ReadNetgen3DMesh(std::istream &input);
    void ReadTrueGridMesh(std::istream &input);
-   void CreateVTKMesh(const Vector &points, const Array<int> &cell_data,
-                      const Array<int> &cell_offsets,
-                      const Array<int> &cell_types,
-                      const Array<int> &cell_attributes,
-                      int &curved, int &read_gf, bool &finalize_topo);
    void ReadVTKMesh(std::istream &input, int &curved, int &read_gf,
                     bool &finalize_topo);
-   void ReadXML_VTKMesh(std::istream &input, int &curved, int &read_gf,
-                        bool &finalize_topo);
    void ReadNURBSMesh(std::istream &input, int &curved, int &read_gf);
    void ReadInlineMesh(std::istream &input, bool generate_edges = false);
    void ReadGmshMesh(std::istream &input, int &curved, int &read_gf);

mesh/mesh_readers.cpp

@@ -12,13 +12,10 @@
 #include "mesh_headers.hpp"
 #include "../fem/fem.hpp"
 #include "../general/text.hpp"
-#include "../general/tinyxml2.h"
 #include "gmsh.hpp"
 
 #include <iostream>
 #include <cstdio>
-#include <vector>
-#include <algorithm>
 
 #ifdef MFEM_USE_NETCDF
 #include "netcdf.h"
@@ -376,398 +373,6 @@ const int Mesh::vtk_quadratic_hex[27] =
    24, 22, 21, 23, 20, 25, 26
 };
 
-void Mesh::CreateVTKMesh(const Vector &points, const Array<int> &cell_data,
-                         const Array<int> &cell_offsets,
-                         const Array<int> &cell_types,
-                         const Array<int> &cell_attributes,
-                         int &curved, int &read_gf, bool &finalize_topo)
-{
-   int np = points.Size()/3;
-   Dim = -1;
-   NumOfElements = cell_types.Size();
-   elements.SetSize(NumOfElements);
-
-   int order = -1;
-   bool legacy_elem = false, lagrange_elem = false;
-
-   int j = 0;
-   for (int i = 0; i < NumOfElements; i++)
-   {
-      int ct = cell_types[i];
-      Geometry::Type geom = VTKGeometry::GetMFEMGeometry(ct);
-      elements[i] = NewElement(geom);
-      if (cell_attributes.Size() > 0)
-      {
-         elements[i]->SetAttribute(cell_attributes[i]);
-      }
-      // VTK ordering of vertices is the same as MFEM ordering of vertices
-      // for all element types *except* prisms, which require a permutation
-      if (geom == Geometry::PRISM && ct != VTKGeometry::LAGRANGE_PRISM)
-      {
-         int prism_vertices[6];
-         for (int k=0; k<6; ++k)
-         {
-            prism_vertices[k] = cell_data[j+VTKGeometry::PrismMap[k]];
-         }
-         elements[i]->SetVertices(prism_vertices);
-      }
-      else
-      {
-         elements[i]->SetVertices(&cell_data[j]);
-      }
-
-      int elem_dim = Geometry::Dimension[geom];
-      int elem_order = VTKGeometry::GetOrder(ct, cell_offsets[i] - j);
-
-      if (VTKGeometry::IsLagrange(ct)) { lagrange_elem = true; }
-      else { legacy_elem = true; }
-
-      MFEM_VERIFY(Dim == -1 || Dim == elem_dim,
-                  "Elements with different dimensions are not supported");
-      MFEM_VERIFY(order == -1 || order == elem_order,
-                  "Elements with different orders are not supported");
-      MFEM_VERIFY(legacy_elem != lagrange_elem,
-                  "Mixing of legacy and Lagrange cell types is not supported");
-      Dim = elem_dim;
-      order = elem_order;
-      j = cell_offsets[i];
-   }
-
-   if (order == 1 && !lagrange_elem)
-   {
-      NumOfVertices = np;
-      vertices.SetSize(np);
-      for (int i = 0; i < np; i++)
-      {
-         vertices[i](0) = points(3*i+0);
-         vertices[i](1) = points(3*i+1);
-         vertices[i](2) = points(3*i+2);
-      }
-      // No boundary is defined in a VTK mesh
-      NumOfBdrElements = 0;
-      FinalizeTopology();
-      CheckElementOrientation(true);
-   }
-   else
-   {
-      // The following section of code is shared for legacy quadratic and the
-      // Lagrange high order elements
-      curved = 1;
-
-      // generate new enumeration for the vertices
-      Array<int> pts_dof(np);
-      pts_dof = -1;
-      // mark vertex points
-      for (int i = 0; i < NumOfElements; i++)
-      {
-         int *v = elements[i]->GetVertices();
-         int nv = elements[i]->GetNVertices();
-         for (int j = 0; j < nv; j++)
-         {
-            if (pts_dof[v[j]] == -1) { pts_dof[v[j]] = 0; }
-         }
-      }
-
-      // The following loop reorders pts_dofs so vertices are visited in
-      // canonical order
-
-      // Keep the original ordering of the vertices
-      int i, n;
-      for (n = i = 0; i < np; i++)
-      {
-         if (pts_dof[i] != -1)
-         {
-            pts_dof[i] = n++;
-         }
-      }
-      // update the element vertices
-      for (int i = 0; i < NumOfElements; i++)
-      {
-         int *v = elements[i]->GetVertices();
-         int nv = elements[i]->GetNVertices();
-         for (int j = 0; j < nv; j++)
-         {
-            v[j] = pts_dof[v[j]];
-         }
-      }
-      // Define the 'vertices' from the 'points' through the 'pts_dof' map
-      NumOfVertices = n;
-      vertices.SetSize(n);
-      for (int i = 0; i < np; i++)
-      {
-         int j = pts_dof[i];
-         if (j != -1)
-         {
-            vertices[j](0) = points(3*i+0);
-            vertices[j](1) = points(3*i+1);
-            vertices[j](2) = points(3*i+2);
-         }
-      }
-
-      // No boundary is defined in a VTK mesh
-      NumOfBdrElements = 0;
-
-      // determine spaceDim based on min/max differences detected each dimension
-      if (vertices.Size() > 0)
-      {
-         double min_value, max_value;
-         for (int d=0; d<3; ++d)
-         {
-            min_value = max_value = vertices[0](d);
-            for (int i = 1; i < vertices.Size(); i++)
-            {
-               min_value = std::min(min_value,vertices[i](d));
-               max_value = std::max(max_value,vertices[i](d));
-               if (min_value != max_value)
-               {
-                  spaceDim++;
-                  break;
-               }
-            }
-         }
-      }
-      // Generate faces and edges so that we can define
-      // FE space on the mesh
-      FinalizeTopology();
-
-      FiniteElementCollection *fec;
-      FiniteElementSpace *fes;
-      if (legacy_elem)
-      {
-         // Define quadratic FE space
-         fec = new QuadraticFECollection;
-         fes = new FiniteElementSpace(this, fec, spaceDim);
-         Nodes = new GridFunction(fes);
-         Nodes->MakeOwner(fec); // Nodes will destroy 'fec' and 'fes'
-         own_nodes = 1;
-
-         // Map vtk points to edge/face/element dofs
-         Array<int> dofs;
-         for (int i = 0; i < NumOfElements; i++)
-         {
-            fes->GetElementDofs(i, dofs);
-            const int *vtk_mfem;
-            switch (elements[i]->GetGeometryType())
-            {
-               case Geometry::TRIANGLE:
-               case Geometry::SQUARE:
-                  vtk_mfem = vtk_quadratic_hex; break; // identity map
-               case Geometry::TETRAHEDRON:
-                  vtk_mfem = vtk_quadratic_tet; break;
-               case Geometry::CUBE:
-                  vtk_mfem = vtk_quadratic_hex; break;
-               case Geometry::PRISM:
-                  vtk_mfem = vtk_quadratic_wedge; break;
-               default:
-                  vtk_mfem = NULL; // suppress a warning
-                  break;
-            }
-
-            int offset = (i == 0) ? 0 : cell_offsets[i-1];
-            for (int j = 0; j < dofs.Size(); j++)
-            {
-               if (pts_dof[cell_data[offset+j]] == -1)
-               {
-                  pts_dof[cell_data[offset+j]] = dofs[vtk_mfem[j]];
-               }
-               else
-               {
-                  if (pts_dof[cell_data[offset+j]] != dofs[vtk_mfem[j]])
-                  {
-                     MFEM_ABORT("VTK mesh: inconsistent quadratic mesh!");
-                  }
-               }
-            }
-         }
-      }
-      else
-      {
-         // Define H1 FE space
-         fec = new H1_FECollection(order,Dim,BasisType::ClosedUniform);
-         fes = new FiniteElementSpace(this, fec, spaceDim);
-         Nodes = new GridFunction(fes);
-         Nodes->MakeOwner(fec); // Nodes will destroy 'fec' and 'fes'
-         own_nodes = 1;
-         Array<int> dofs;
-
-         std::map<Geometry::Type,Array<int>> vtk_inv_maps;
-         std::map<Geometry::Type,const Array<int>*> lex_orderings;
-
-         int i, n;
-         for (n = i = 0; i < NumOfElements; i++)
-         {
-            Geometry::Type geom = GetElementBaseGeometry(i);
-            fes->GetElementDofs(i, dofs);
-
-            Array<int> &vtk_inv_map = vtk_inv_maps[geom];
-            if (vtk_inv_map.Size() == 0)
-            {
-               Array<int> vtk_map;
-               CreateVTKElementConnectivity(vtk_map, geom, order);
-               vtk_inv_map.SetSize(vtk_map.Size());
-               for (int j=0; j<vtk_map.Size(); ++j)
-               {
-                  vtk_inv_map[vtk_map[j]] = j;
-               }
-            }
-            const Array<int> *&lex_ordering = lex_orderings[geom];
-            if (!lex_ordering)
-            {
-               const FiniteElement *fe = fes->GetFE(i);
-               const NodalFiniteElement *nodal_fe =
-                  dynamic_cast<const NodalFiniteElement*>(fe);
-               MFEM_ASSERT(nodal_fe != NULL, "Unsupported element type");
-               lex_ordering = &nodal_fe->GetLexicographicOrdering();
-            }
-
-            for (int lex_idx = 0; lex_idx < dofs.Size(); lex_idx++)
-            {
-               int mfem_idx = (*lex_ordering)[lex_idx];
-               int vtk_idx = vtk_inv_map[lex_idx];
-               int pt_idx = cell_data[n + vtk_idx];
-               if (pts_dof[pt_idx] == -1)
-               {
-                  pts_dof[pt_idx] = dofs[mfem_idx];
-               }
-               else
-               {
-                  if (pts_dof[pt_idx] != dofs[mfem_idx])
-                  {
-                     MFEM_ABORT("VTK mesh: inconsistent Lagrange mesh!");
-                  }
-               }
-            }
-            n += dofs.Size();
-         }
-      }
-      // Define the 'Nodes' from the 'points' through the 'pts_dof' map
-      Array<int> dofs;
-      for (int i = 0; i < np; i++)
-      {
-         dofs.SetSize(1);
-         if (pts_dof[i] != -1)
-         {
-            dofs[0] = pts_dof[i];
-            fes->DofsToVDofs(dofs);
-            for (int d = 0; d < dofs.Size(); d++)
-            {
-               (*Nodes)(dofs[d]) = points(3*i+d);
-            }
-         }
-      }
-      read_gf = 0;
-   }
-}
-
-void Mesh::ReadXML_VTKMesh(std::istream &input, int &curved, int &read_gf,
-                           bool &finalize_topo)
-{
-   using namespace tinyxml2;
-
-   const char *erstr = "XML parsing error";
-
-   // Read entire stream into buffer
-   std::istreambuf_iterator<char> eos;
-   std::vector<char> buf(std::istreambuf_iterator<char>(input), eos);
-   buf.push_back('\0'); // null-terminate buffer
-
-   XMLDocument xml;
-   xml.Parse(buf.data());
-   MFEM_VERIFY(xml.ErrorID() == XML_SUCCESS, erstr);
-
-   const XMLElement *vtkfile = xml.FirstChildElement();
-   MFEM_VERIFY(vtkfile, erstr);
-   MFEM_VERIFY(std::string(vtkfile->Name()) == "VTKFile", erstr);
-   const XMLElement *vtu = vtkfile->FirstChildElement();
-   MFEM_VERIFY(vtu, erstr);
-   MFEM_VERIFY(std::string(vtu->Name()) == "UnstructuredGrid", erstr);
-
-   // Count the number of points and cells
-   const XMLElement *piece = vtu->FirstChildElement();
-   MFEM_VERIFY(std::string(piece->Name()) == "Piece", erstr);
-   MFEM_VERIFY(piece->NextSiblingElement() == NULL,
-               "XML VTK meshes with more than one Piece are not supported");
-   int npts = piece->IntAttribute("NumberOfPoints");
-   int ncells = piece->IntAttribute("NumberOfCells");
-
-   // Read the points
-   Vector points(3*npts);
-   const XMLElement *pts_xml;
-   for (pts_xml = piece->FirstChildElement();
-        pts_xml != NULL;
-        pts_xml = pts_xml->NextSiblingElement())
-   {
-      if (std::string(pts_xml->Name()) == "Points")
-      {
-         const XMLElement *pts_data = pts_xml->FirstChildElement();
-         MFEM_VERIFY(std::string(pts_data->Name()) == "DataArray", erstr);
-         MFEM_VERIFY(std::string(pts_data->Attribute("Name")) == "Points",
-                     erstr);
-         MFEM_VERIFY(pts_data->IntAttribute("NumberOfComponents") == 3,
-                     "XML VTK Points DataArray must have 3 components");
-         const char *pts_txt = pts_data->GetText();
-         MFEM_VERIFY(pts_txt != NULL, erstr);
-         std::istringstream pts_stream(pts_txt);
-         points.Load(pts_stream, 3*npts);
-         break;
-      }
-   }
-   if (pts_xml == NULL) { MFEM_ABORT(erstr); }
-
-   // Read the cells
-   Array<int> cell_data, cell_offsets, cell_types;
-   const XMLElement *cells_xml;
-   for (cells_xml = piece->FirstChildElement();
-        cells_xml != NULL;
-        cells_xml = cells_xml->NextSiblingElement())
-   {
-      if (std::string(cells_xml->Name()) == "Cells")
-      {
-         const char *cell_data_txt = NULL;
-         for (const XMLElement *data_xml = cells_xml->FirstChildElement();
-              data_xml != NULL;
-              data_xml = data_xml->NextSiblingElement())
-         {
-            MFEM_VERIFY(std::string(data_xml->Name()) == "DataArray", erstr);
-            std::string data_name(data_xml->Attribute("Name"));
-            const char *data_txt = data_xml->GetText();
-            MFEM_VERIFY(data_txt != NULL, erstr);
-            if (data_name == "offsets")
-            {
-               std::istringstream data_stream(data_txt);
-               cell_offsets.Load(ncells, data_stream);
-            }
-            else if (data_name == "types")
-            {
-               std::istringstream data_stream(data_txt);
-               cell_types.Load(ncells, data_stream);
-            }
-            else if (data_name == "connectivity")
-            {
-               // Have to read the connectivity after the offsets, because we
-               // don't know how many points to read until we have the offsets
-               // (size of connectivity array is equal to the last offset), so
-               // store the data pointer and read this array later.
-               cell_data_txt = data_txt;
-            }
-         }
-         MFEM_VERIFY(cell_offsets.Size() == ncells, erstr);
-         MFEM_VERIFY(cell_types.Size() == ncells, erstr);
-         MFEM_VERIFY(cell_data_txt != NULL, erstr);
-         int cell_data_size = cell_offsets.Last();
-         std::istringstream cell_data_stream(cell_data_txt);
-         cell_data.Load(cell_data_size, cell_data_stream);
-         break;
-      }
-   }
-   if (cells_xml == NULL) { MFEM_ABORT(erstr); }
-
-   // Currently don't support reading cell attributes from VTK mesh
-   Array<int> cell_attributes;
-   CreateVTKMesh(points, cell_data, cell_offsets, cell_types, cell_attributes,
-                 curved, read_gf, finalize_topo);
-}
-
 void Mesh::ReadVTKMesh(std::istream &input, int &curved, int &read_gf,
                        bool &finalize_topo)
 {
@@ -777,6 +382,8 @@ void Mesh::ReadVTKMesh(std::istream &input, int &curved, int &read_gf,
    //   * https://lorensen.github.io/VTKExamples/site/VTKFileFormats
    //   * https://www.kitware.com/products/books/VTKUsersGuide.pdf
 
+   int i, j, n, attr;
+
    string buff;
    getline(input, buff); // comment line
    getline(input, buff);
@@ -805,78 +412,141 @@ void Mesh::ReadVTKMesh(std::istream &input, int &curved, int &read_gf,
       }
    }
    while (buff != "POINTS");
-
+   int np = 0;
    Vector points;
-   int np;
-   input >> np >> ws;
-   getline(input, buff); // "double"
-   points.Load(input, 3*np);
-
-   //skip metadata
-   // Looks like:
-   // METADATA
-   //INFORMATION 2
-   //NAME L2_NORM_RANGE LOCATION vtkDataArray
-   //DATA 2 0 5.19615
-   //NAME L2_NORM_FINITE_RANGE LOCATION vtkDataArray
-   //DATA 2 0 5.19615
-   do
    {
-      input >> buff;
-      if (!input.good())
+      input >> np >> ws;
+      points.SetSize(3*np);
+      getline(input, buff); // "double"
+      for (i = 0; i < points.Size(); i++)
       {
-         MFEM_ABORT("VTK mesh does not have CELLS data!");
+         input >> points(i);
       }
    }
-   while (buff != "CELLS");
 
    // Read the cells
-   Array<int> cell_data, cell_offsets;
+   NumOfElements = n = 0;
+   Array<int> cells_data;
+   input >> ws >> buff;
    if (buff == "CELLS")
    {
-      int ncells, n;
-      input >> ncells >> n >> ws;
-      cell_offsets.SetSize(ncells);
-      cell_data.SetSize(n - ncells);
-      int offset = 0;
-      for (int i=0; i<ncells; ++i)
+      input >> NumOfElements >> n >> ws;
+      cells_data.SetSize(n);
+      for (i = 0; i < n; i++)
       {
-         int nv;
-         input >> nv;
-         cell_offsets[i] = offset + nv;
-         for (int j=0; j<nv; ++j)
-         {
-            input >> cell_data[offset + j];
-         }
-         offset += nv;
+         input >> cells_data[i];
       }
    }
 
    // Read the cell types
+   Dim = -1;
+   int order = -1;
    input >> ws >> buff;
-   Array<int> cell_types;
-   int ncells;
    if (buff == "CELL_TYPES")
    {
-      input >> ncells;
-      cell_types.Load(ncells, input);
+      input >> NumOfElements;
+      elements.SetSize(NumOfElements);
+      for (j = i = 0; i < NumOfElements; i++)
+      {
+         int ct, elem_dim, elem_order = 1;
+         input >> ct;
+         switch (ct)
+         {
+            case 5:   // triangle
+               elem_dim = 2;
+               elements[i] = new Triangle(&cells_data[j+1]);
+               break;
+            case 9:   // quadrilateral
+               elem_dim = 2;
+               elements[i] = new Quadrilateral(&cells_data[j+1]);
+               break;
+            case 10:  // tetrahedron
+               elem_dim = 3;
+#ifdef MFEM_USE_MEMALLOC
+               elements[i] = TetMemory.Alloc();
+               elements[i]->SetVertices(&cells_data[j+1]);
+#else
+               elements[i] = new Tetrahedron(&cells_data[j+1]);
+#endif
+               break;
+            case 12:  // hexahedron
+               elem_dim = 3;
+               elements[i] = new Hexahedron(&cells_data[j+1]);
+               break;
+            case 13:  // wedge
+               elem_dim = 3;
+               // switch between vtk vertex ordering and mfem vertex ordering:
+               // swap vertices (1,2) and (4,5)
+               elements[i] =
+                  new Wedge(cells_data[j+1], cells_data[j+3], cells_data[j+2],
+                            cells_data[j+4], cells_data[j+6], cells_data[j+5]);
+               break;
+
+            case 22:  // quadratic triangle
+               elem_dim = 2;
+               elem_order = 2;
+               elements[i] = new Triangle(&cells_data[j+1]);
+               break;
+            case 28:  // biquadratic quadrilateral
+               elem_dim = 2;
+               elem_order = 2;
+               elements[i] = new Quadrilateral(&cells_data[j+1]);
+               break;
+            case 24:  // quadratic tetrahedron
+               elem_dim = 3;
+               elem_order = 2;
+#ifdef MFEM_USE_MEMALLOC
+               elements[i] = TetMemory.Alloc();
+               elements[i]->SetVertices(&cells_data[j+1]);
+#else
+               elements[i] = new Tetrahedron(&cells_data[j+1]);
+#endif
+               break;
+            case 32: // biquadratic-quadratic wedge
+               elem_dim = 3;
+               elem_order = 2;
+               // switch between vtk vertex ordering and mfem vertex ordering:
+               // swap vertices (1,2) and (4,5)
+               elements[i] =
+                  new Wedge(cells_data[j+1], cells_data[j+3], cells_data[j+2],
+                            cells_data[j+4], cells_data[j+6], cells_data[j+5]);
+               break;
+            case 29:  // triquadratic hexahedron
+               elem_dim = 3;
+               elem_order = 2;
+               elements[i] = new Hexahedron(&cells_data[j+1]);
+               break;
+            default:
+               MFEM_ABORT("VTK mesh : cell type " << ct << " is not supported!");
+               return;
+         }
+         MFEM_VERIFY(Dim == -1 || Dim == elem_dim,
+                     "elements with different dimensions are not supported");
+         MFEM_VERIFY(order == -1 || order == elem_order,
+                     "elements with different orders are not supported");
+         Dim = elem_dim;
+         order = elem_order;
+         j += cells_data[j] + 1;
+      }
    }
 
-   // Read cell attributes
+   // Read attributes
    streampos sp = input.tellg();
    input >> ws >> buff;
-   Array<int> cell_attributes;
    if (buff == "CELL_DATA")
    {
-      int n;
       input >> n >> ws;
       getline(input, buff);
       filter_dos(buff);
       // "SCALARS material dataType numComp"
-      if (buff.rfind("SCALARS material") == 0)
+      if (!strncmp(buff.c_str(), "SCALARS material", 16))
       {
          getline(input, buff); // "LOOKUP_TABLE default"
-         cell_attributes.Load(ncells, input);
+         for (i = 0; i < NumOfElements; i++)
+         {
+            input >> attr;
+            elements[i]->SetAttribute(attr);
+         }
       }
       else
       {
@@ -888,9 +558,138 @@ void Mesh::ReadVTKMesh(std::istream &input, int &curved, int &read_gf,
       input.seekg(sp);
    }
 
-   CreateVTKMesh(points, cell_data, cell_offsets, cell_types, cell_attributes,
-                 curved, read_gf, finalize_topo);
-} // end ReadVTKMesh
+   if (order == 1)
+   {
+      cells_data.DeleteAll();
+      NumOfVertices = np;
+      vertices.SetSize(np);
+      for (i = 0; i < np; i++)
+      {
+         vertices[i](0) = points(3*i+0);
+         vertices[i](1) = points(3*i+1);
+         vertices[i](2) = points(3*i+2);
+      }
+      points.Destroy();
+
+      // No boundary is defined in a VTK mesh
+      NumOfBdrElements = 0;
+   }
+   else if (order == 2)
+   {
+      curved = 1;
+
+      // generate new enumeration for the vertices
+      Array<int> pts_dof(np);
+      pts_dof = -1;
+      for (n = i = 0; i < NumOfElements; i++)
+      {
+         int *v = elements[i]->GetVertices();
+         int nv = elements[i]->GetNVertices();
+         for (j = 0; j < nv; j++)
+            if (pts_dof[v[j]] == -1)
+            {
+               pts_dof[v[j]] = n++;
+            }
+      }
+      // keep the original ordering of the vertices
+      for (n = i = 0; i < np; i++)
+         if (pts_dof[i] != -1)
+         {
+            pts_dof[i] = n++;
+         }
+      // update the element vertices
+      for (i = 0; i < NumOfElements; i++)
+      {
+         int *v = elements[i]->GetVertices();
+         int nv = elements[i]->GetNVertices();
+         for (j = 0; j < nv; j++)
+         {
+            v[j] = pts_dof[v[j]];
+         }
+      }
+      // Define the 'vertices' from the 'points' through the 'pts_dof' map
+      NumOfVertices = n;
+      vertices.SetSize(n);
+      for (i = 0; i < np; i++)
+      {
+         if ((j = pts_dof[i]) != -1)
+         {
+            vertices[j](0) = points(3*i+0);
+            vertices[j](1) = points(3*i+1);
+            vertices[j](2) = points(3*i+2);
+         }
+      }
+
+      // No boundary is defined in a VTK mesh
+      NumOfBdrElements = 0;
+
+      // Generate faces and edges so that we can define quadratic
+      // FE space on the mesh
+      FinalizeTopology();
+      finalize_topo = false;
+
+      // Define quadratic FE space
+      FiniteElementCollection *fec = new QuadraticFECollection;
+      FiniteElementSpace *fes = new FiniteElementSpace(this, fec, Dim);
+      Nodes = new GridFunction(fes);
+      Nodes->MakeOwner(fec); // Nodes will destroy 'fec' and 'fes'
+      own_nodes = 1;
+
+      // Map vtk points to edge/face/element dofs
+      Array<int> dofs;
+      for (n = i = 0; i < NumOfElements; i++)
+      {
+         fes->GetElementDofs(i, dofs);
+         const int *vtk_mfem;
+         switch (elements[i]->GetGeometryType())
+         {
+            case Geometry::TRIANGLE:
+            case Geometry::SQUARE:
+               vtk_mfem = vtk_quadratic_hex; break; // identity map
+            case Geometry::TETRAHEDRON:
+               vtk_mfem = vtk_quadratic_tet; break;
+            case Geometry::CUBE:
+               vtk_mfem = vtk_quadratic_hex; break;
+            case Geometry::PRISM:
+               vtk_mfem = vtk_quadratic_wedge; break;
+            default:
+               vtk_mfem = NULL; // suppress a warning
+               break;
+         }
+
+         for (n++, j = 0; j < dofs.Size(); j++, n++)
+         {
+            if (pts_dof[cells_data[n]] == -1)
+            {
+               pts_dof[cells_data[n]] = dofs[vtk_mfem[j]];
+            }
+            else
+            {
+               if (pts_dof[cells_data[n]] != dofs[vtk_mfem[j]])
+               {
+                  MFEM_ABORT("VTK mesh : inconsistent quadratic mesh!");
+               }
+            }
+         }
+      }
+
+      // Define the 'Nodes' from the 'points' through the 'pts_dof' map
+      for (i = 0; i < np; i++)
+      {
+         dofs.SetSize(1);
+         if ((dofs[0] = pts_dof[i]) != -1)
+         {
+            fes->DofsToVDofs(dofs);
+            for (j = 0; j < dofs.Size(); j++)
+            {
+               (*Nodes)(dofs[j]) = points(3*i+j);
+            }
+         }
+      }
+
+      read_gf = 0;
+   }
+}
 
 void Mesh::ReadNURBSMesh(std::istream &input, int &curved, int &read_gf)
 {
@@ -916,8 +715,8 @@ void Mesh::ReadNURBSMesh(std::istream &input, int &curved, int &read_gf)
       NURBSext->SetCoordsFromPatches(*Nodes);
       own_nodes = 1;
       read_gf = 0;
-      spaceDim = Nodes->VectorDim();
-      for (int i = 0; i < spaceDim; i++)
+      int vd = Nodes->VectorDim();
+      for (int i = 0; i < vd; i++)
       {
          Vector vert_val;
          Nodes->GetNodalValues(vert_val, i+1);

mesh/vtk.cpp

@@ -18,137 +18,7 @@
 namespace mfem
 {
 
-const int VTKGeometry::Map[Geometry::NUM_GEOMETRIES] =
-{
-   POINT, SEGMENT, TRIANGLE, SQUARE, TETRAHEDRON, CUBE, PRISM
-};
-
-const int VTKGeometry::QuadraticMap[Geometry::NUM_GEOMETRIES] =
-{
-   POINT, QUADRATIC_SEGMENT, QUADRATIC_TRIANGLE, BIQUADRATIC_SQUARE,
-   QUADRATIC_TETRAHEDRON, TRIQUADRATIC_CUBE, BIQUADRATIC_QUADRATIC_PRISM
-};
-
-const int VTKGeometry::HighOrderMap[Geometry::NUM_GEOMETRIES] =
-{
-   POINT, LAGRANGE_SEGMENT, LAGRANGE_TRIANGLE, LAGRANGE_SQUARE,
-   LAGRANGE_TETRAHEDRON, LAGRANGE_CUBE, LAGRANGE_PRISM
-};
-
-const int VTKGeometry::PrismMap[6] = {0, 2, 1, 3, 5, 4};
-
-const int *VTKGeometry::VertexPermutation[Geometry::NUM_GEOMETRIES] =
-{
-   NULL, NULL, NULL, NULL, NULL, NULL, VTKGeometry::PrismMap
-};
-
-Geometry::Type VTKGeometry::GetMFEMGeometry(int vtk_geom)
-{
-   switch (vtk_geom)
-   {
-      case POINT:
-         return Geometry::POINT;
-      case SEGMENT:
-      case QUADRATIC_SEGMENT:
-      case LAGRANGE_SEGMENT:
-         return Geometry::SEGMENT;
-      case TRIANGLE:
-      case QUADRATIC_TRIANGLE:
-      case LAGRANGE_TRIANGLE:
-         return Geometry::TRIANGLE;
-      case SQUARE:
-      case BIQUADRATIC_SQUARE:
-      case LAGRANGE_SQUARE:
-         return Geometry::SQUARE;
-      case TETRAHEDRON:
-      case QUADRATIC_TETRAHEDRON:
-      case LAGRANGE_TETRAHEDRON:
-         return Geometry::TETRAHEDRON;
-      case CUBE:
-      case TRIQUADRATIC_CUBE:
-      case LAGRANGE_CUBE:
-         return Geometry::CUBE;
-      case PRISM:
-      case BIQUADRATIC_QUADRATIC_PRISM:
-      case LAGRANGE_PRISM:
-         return Geometry::PRISM;
-      default:
-         return Geometry::INVALID;
-   }
-}
-
-bool VTKGeometry::IsLagrange(int vtk_geom)
-{
-   return vtk_geom >= LAGRANGE_SEGMENT && vtk_geom <= LAGRANGE_PRISM;
-}
-
-bool VTKGeometry::IsQuadratic(int vtk_geom)
-{
-   return vtk_geom >= QUADRATIC_SEGMENT
-          && vtk_geom <= BIQUADRATIC_QUADRATIC_PRISM;
-}
-
-int VTKGeometry::GetOrder(int vtk_geom, int npoints)
-{
-   if (IsQuadratic(vtk_geom))
-   {
-      return 2;
-   }
-   else if (IsLagrange(vtk_geom))
-   {
-      switch (vtk_geom)
-      {
-         case LAGRANGE_SEGMENT:
-            return npoints - 1;
-         case LAGRANGE_TRIANGLE:
-            return (std::sqrt(8*npoints + 1) - 3)/2;
-         case LAGRANGE_SQUARE:
-            return std::round(std::sqrt(npoints)) - 1;
-         case LAGRANGE_TETRAHEDRON:
-            switch (npoints)
-            {
-               // Note that for given order, npoints is given by
-               // npoints_order = (order + 1)*(order + 2)*(order + 3)/6,
-               case 4: return 1;
-               case 10: return 2;
-               case 20: return 3;
-               case 35: return 4;
-               case 56: return 5;
-               case 84: return 6;
-               case 120: return 7;
-               case 165: return 8;
-               case 220: return 9;
-               case 286: return 10;
-               default:
-               {
-                  constexpr int max_order = 20;
-                  int order = 11, npoints_order;
-                  for (; order<max_order; ++order)
-                  {
-                     npoints_order = (order + 1)*(order + 2)*(order + 3)/6;
-                     if (npoints_order == npoints) { break; }
-                  }
-                  MFEM_VERIFY(npoints == npoints_order, "");
-                  return order;
-               }
-            }
-         case LAGRANGE_CUBE:
-            return std::round(std::cbrt(npoints)) - 1;
-         case LAGRANGE_PRISM:
-         {
-            const double n = npoints;
-            static const double third = 1.0/3.0;
-            static const double ninth = 1.0/9.0;
-            static const double twentyseventh = 1.0/27.0;
-            const double term =
-               std::cbrt(third*sqrt(third)*sqrt((27.0*n - 2.0)*n) + n
-                         - twentyseventh);
-            return std::round(term + ninth / term - 4*third);
-         }
-      }
-   }
-   return 1;
-}
+const int vtk_prism_perm[6] = {0, 2, 1, 3, 5, 4};
 
 int BarycentricToVTKTriangle(int *b, int ref)
 {
@@ -509,12 +379,15 @@ void CreateVTKElementConnectivity(Array<int> &con, Geometry::Type geom, int ref)
    {
       int idx = 0;
       int b[4];
-      for (b[2]=0; b[2]<=ref; b[2]++)
+      for (int k=0; k<=ref; k++)
       {
-         for (b[1]=0; b[1]<=ref-b[2]; b[1]++)
+         for (int j=0; j<=k; j++)
          {
-            for (b[0]=0; b[0]<=ref-b[1]-b[2]; b[0]++)
+            for (int i=0; i<=j; i++)
             {
+               b[0] = k-j;
+               b[1] = i;
+               b[2] = j-i;
                b[3] = ref-b[0]-b[1]-b[2];
                con[BarycentricToVTKTetra(b, ref)] = idx++;
             }

mesh/vtk.hpp

@@ -17,46 +17,7 @@
 namespace mfem
 {
 
-// Helpers for reading and writing VTK format
-
-// VTK element types defined at: https://git.io/JvZLm
-struct VTKGeometry
-{
-   static const int POINT = 1;
-   static const int SEGMENT = 3;
-   static const int TRIANGLE = 5;
-   static const int SQUARE = 9;
-   static const int TETRAHEDRON = 10;
-   static const int CUBE = 12;
-   static const int PRISM = 13;
-
-   static const int QUADRATIC_SEGMENT = 21;
-   static const int QUADRATIC_TRIANGLE = 22;
-   static const int BIQUADRATIC_SQUARE = 28;
-   static const int QUADRATIC_TETRAHEDRON = 24;
-   static const int TRIQUADRATIC_CUBE = 29;
-   static const int QUADRATIC_PRISM = 26;
-   static const int BIQUADRATIC_QUADRATIC_PRISM = 32;
-
-   static const int LAGRANGE_SEGMENT = 68;
-   static const int LAGRANGE_TRIANGLE = 69;
-   static const int LAGRANGE_SQUARE = 70;
-   static const int LAGRANGE_TETRAHEDRON = 71;
-   static const int LAGRANGE_CUBE = 72;
-   static const int LAGRANGE_PRISM = 73;
-
-   static const int PrismMap[6];
-   static const int *VertexPermutation[Geometry::NUM_GEOMETRIES];
-
-   static const int Map[Geometry::NUM_GEOMETRIES];
-   static const int QuadraticMap[Geometry::NUM_GEOMETRIES];
-   static const int HighOrderMap[Geometry::NUM_GEOMETRIES];
-
-   static Geometry::Type GetMFEMGeometry(int vtk_geom);
-   static bool IsLagrange(int vtk_geom);
-   static bool IsQuadratic(int vtk_geom);
-   static int GetOrder(int vtk_geom, int npoints);
-};
+// Helpers for writing to the VTK format
 
 enum class VTKFormat
 {
@@ -65,6 +26,8 @@ enum class VTKFormat
    BINARY32
 };
 
+extern const int vtk_prism_perm[6];
+
 /// Create the VTK element connectivity array for a given element geometry and
 /// refinement level. Converts node numbers from MFEM to VTK ordering.
 void CreateVTKElementConnectivity(Array<int> &con, Geometry::Type geom,

tests/unit/CMakeLists.txt

@@ -42,7 +42,6 @@ set(UNIT_TESTS_SRCS
   fem/test_2d_bilininteg.cpp
   fem/test_3d_bilininteg.cpp
   fem/test_assemblediagonalpa.cpp
-  fem/test_blocknonlinearform.cpp
   fem/test_calcshape.cpp
   fem/test_datacollection.cpp
   fem/test_estimator.cpp
@@ -51,7 +50,6 @@ set(UNIT_TESTS_SRCS
   fem/test_intrules.cpp
   fem/test_intruletypes.cpp
   fem/test_inversetransform.cpp
-  fem/test_lexicographic_ordering.cpp
   fem/test_lin_interp.cpp
   fem/test_linear_fes.cpp
   fem/test_operatorjacobismoother.cpp
@@ -61,7 +59,7 @@ set(UNIT_TESTS_SRCS
   fem/test_pa_idinterp.cpp
   fem/test_quadf_coef.cpp
   fem/test_quadraturefunc.cpp
-  fem/test_sum_bilin.cpp
+  fem/test_blocknonlinearform.cpp
   miniapps/test_sedov.cpp
 )
 

tests/unit/fem/test_lexicographic_ordering.cpp

@@ -1,64 +0,0 @@
-// Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
-// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
-// LICENSE and NOTICE for details. LLNL-CODE-806117.
-//
-// This file is part of the MFEM library. For more information and source code
-// availability visit https://mfem.org.
-//
-// MFEM is free software; you can redistribute it and/or modify it under the
-// terms of the BSD-3 license. We welcome feedback and contributions, see file
-// CONTRIBUTING.md for details.
-
-#include "mfem.hpp"
-#include "unit_tests.hpp"
-
-using namespace mfem;
-
-void VerifyOrdering(NodalFiniteElement &el)
-{
-   int order = el.GetOrder();
-   Geometry::Type geom = el.GetGeomType();
-   const Array<int> &p = el.GetLexicographicOrdering();
-
-   GeometryRefiner refiner;
-   refiner.SetType(BasisType::GaussLobatto);
-   RefinedGeometry *ref_geom = refiner.Refine(geom, order);
-
-   double error = 0.0;
-
-   for (int i=0; i<el.GetDof(); ++i)
-   {
-      int pi = (p.Size() > 0) ? p[i] : i;
-      error += std::fabs(el.GetNodes()[pi].x - ref_geom->RefPts[i].x);
-      error += std::fabs(el.GetNodes()[pi].y - ref_geom->RefPts[i].y);
-      error += std::fabs(el.GetNodes()[pi].z - ref_geom->RefPts[i].z);
-   }
-
-   REQUIRE(error == MFEM_Approx(0.0));
-}
-
-template <typename T> void VerifyOrdering(int order)
-{
-   T el(order, BasisType::GaussLobatto);
-   Geometry::Type geom = el.GetGeomType();
-   INFO("order " << order << " " << Geometry::Name[geom]);
-   VerifyOrdering(el);
-}
-
-TEST_CASE("Lexicographic Ordering", "[FiniteElement,Geometry]")
-{
-   auto order = GENERATE(1, 2, 3, 4, 5, 6);
-   VerifyOrdering<H1_SegmentElement>(order);
-   VerifyOrdering<H1_TriangleElement>(order);
-   VerifyOrdering<H1_QuadrilateralElement>(order);
-   VerifyOrdering<H1_TetrahedronElement>(order);
-   VerifyOrdering<H1_HexahedronElement>(order);
-   VerifyOrdering<H1_WedgeElement>(order);
-
-   VerifyOrdering<L2_SegmentElement>(order);
-   VerifyOrdering<L2_TriangleElement>(order);
-   VerifyOrdering<L2_QuadrilateralElement>(order);
-   VerifyOrdering<L2_TetrahedronElement>(order);
-   VerifyOrdering<L2_HexahedronElement>(order);
-   VerifyOrdering<L2_WedgeElement>(order);
-}

tests/unit/fem/test_sum_bilin.cpp

@@ -1,231 +0,0 @@
-// Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
-// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
-// LICENSE and NOTICE for details. LLNL-CODE-806117.
-//
-// This file is part of the MFEM library. For more information and source code
-// availability visit https://mfem.org.
-//
-// MFEM is free software; you can redistribute it and/or modify it under the
-// terms of the BSD-3 license. We welcome feedback and contributions, see file
-// CONTRIBUTING.md for details.
-
-#include "unit_tests.hpp"
-#include "mfem.hpp"
-#include <fstream>
-#include <iostream>
-
-using namespace mfem;
-
-namespace pa_kernels
-{
-
-TEST_CASE("H1 SumIntegrator", "[SumIntegrator][PartialAssembly]")
-{
-   Mesh mesh(1, 1, 1, Element::HEXAHEDRON);
-   H1_FECollection fec(2, mesh.Dimension());
-   FiniteElementSpace fes(&mesh, &fec);
-
-   MassIntegrator integ1;
-   DiffusionIntegrator integ2;
-
-   SumIntegrator integ_sum(true);
-   integ_sum.AddIntegrator(new MassIntegrator);
-   integ_sum.AddIntegrator(new DiffusionIntegrator);
-
-   const FiniteElement &el = *fes.GetFE(0);
-   ElementTransformation &T = *mesh.GetElementTransformation(0);
-   DenseMatrix m1, m_tmp, m2;
-
-   // AssembleElementMatrix
-   integ1.AssembleElementMatrix(el, T, m1);
-   integ2.AssembleElementMatrix(el, T, m_tmp);
-   m1 += m_tmp;
-   integ_sum.AssembleElementMatrix(el, T, m2);
-   m1 -= m2;
-   REQUIRE(m1.MaxMaxNorm() == MFEM_Approx(0.0));
-
-   // AssembleElementMatrix2
-   integ1.AssembleElementMatrix2(el, el, T, m1);
-   integ2.AssembleElementMatrix2(el, el, T, m_tmp);
-   m1 += m_tmp;
-   integ_sum.AssembleElementMatrix2(el, el, T, m2);
-   m1 -= m2;
-   REQUIRE(m1.MaxMaxNorm() == MFEM_Approx(0.0));
-
-   // PA
-   integ1.AssemblePA(fes);
-   integ2.AssemblePA(fes);
-   integ_sum.AssemblePA(fes);
-
-   int n = fes.GetTrueVSize();
-   Vector x(n), y1(n), y2(n);
-   Vector diag1(n), diag_tmp(n), diag2(n);
-   x.Randomize(1);
-
-   // AddMultPA
-   y1 = 0.0;
-   y2 = 0.0;
-   integ1.AddMultPA(x, y1);
-   integ2.AddMultPA(x, y1);
-   integ_sum.AddMultPA(x, y2);
-   y1 -= y2;
-   REQUIRE(y1.Normlinf() == MFEM_Approx(0.0));
-
-   // AddMultTransposePA
-   y1 = 0.0;
-   y2 = 0.0;
-   integ1.AddMultTransposePA(x, y1);
-   integ2.AddMultTransposePA(x, y1);
-   integ_sum.AddMultTransposePA(x, y2);
-   y1 -= y2;
-   REQUIRE(y1.Normlinf() == MFEM_Approx(0.0));
-
-   // AssembleDiagonalPA
-   diag1 = 0.0;
-   diag_tmp = 0.0;
-   diag2 = 0.0;
-   integ1.AssembleDiagonalPA(diag1);
-   integ2.AssembleDiagonalPA(diag_tmp);
-   diag1 += diag_tmp;
-   integ_sum.AssembleDiagonalPA(diag2);
-   diag1 -= diag2;
-   REQUIRE(diag1.Normlinf() == MFEM_Approx(0.0));
-
-   // MF
-#ifdef MFEM_USE_CEED
-   if (DeviceCanUseCeed())
-   {
-      integ1.AssembleMF(fes);
-      integ2.AssembleMF(fes);
-      integ_sum.AssembleMF(fes);
-
-      // AddMultMF
-      y1 = 0.0;
-      y2 = 0.0;
-      integ1.AddMultMF(x, y1);
-      integ2.AddMultMF(x, y1);
-      integ_sum.AddMultMF(x, y2);
-      y1 -= y2;
-      REQUIRE(y1.Normlinf() == MFEM_Approx(0.0));
-
-      // AddMultTransposeMF
-      y1 = 0.0;
-      y2 = 0.0;
-      integ1.AddMultTransposeMF(x, y1);
-      integ2.AddMultTransposeMF(x, y1);
-      integ_sum.AddMultTransposeMF(x, y2);
-      y1 -= y2;
-      REQUIRE(y1.Normlinf() == MFEM_Approx(0.0));
-
-      // AssembleDiagonalMF
-      integ1.AssembleDiagonalMF(diag1);
-      integ2.AssembleDiagonalMF(diag_tmp);
-      diag1 += diag_tmp;
-      integ_sum.AssembleDiagonalMF(diag2);
-      diag1 -= diag2;
-      REQUIRE(diag1.Normlinf() == MFEM_Approx(0.0));
-   }
-#endif
-}
-
-TEST_CASE("DG SumIntegrator", "[SumIntegrator][PartialAssembly]")
-{
-   Mesh mesh(2, 1, 1, Element::HEXAHEDRON);
-   DG_FECollection fec(2, mesh.Dimension(), BasisType::GaussLobatto);
-   FiniteElementSpace fes(&mesh, &fec);
-
-   Vector v(mesh.Dimension());
-   v = 1.0;
-   VectorConstantCoefficient v_coeff(v);
-
-   DGTraceIntegrator integ1(v_coeff, 1.0, 2.0);
-   DGTraceIntegrator integ2(v_coeff, 3.0, 4.0);
-
-   SumIntegrator integ_sum(true);
-   integ_sum.AddIntegrator(new DGTraceIntegrator(v_coeff, 1.0, 2.0));
-   integ_sum.AddIntegrator(new DGTraceIntegrator(v_coeff, 3.0, 4.0));
-
-   DenseMatrix m1, m_tmp, m2;
-
-   // AssembleFaceMatrix
-   int nfaces = mesh.GetNumFaces();
-   for (int i = 0; i < nfaces; i++)
-   {
-      FaceElementTransformations *tr = mesh.GetFaceElementTransformations(i);
-      const FiniteElement &el0 = *fes.GetFE(tr->Elem1No);
-      const FiniteElement &el1 = (tr->Elem2No >= 0) ? *fes.GetFE(tr->Elem2No) : el0;
-      integ1.AssembleFaceMatrix(el0, el1, *tr, m1);
-      integ2.AssembleFaceMatrix(el0, el1, *tr, m_tmp);
-      m1 += m_tmp;
-      integ_sum.AssembleFaceMatrix(el0, el1, *tr, m2);
-      m1 -= m2;
-      REQUIRE(m1.MaxMaxNorm() == MFEM_Approx(0.0));
-   }
-
-   // PA interior
-   integ1.AssemblePAInteriorFaces(fes);
-   integ2.AssemblePAInteriorFaces(fes);
-   integ_sum.AssemblePAInteriorFaces(fes);
-
-   const Operator *R_int = fes.GetFaceRestriction(
-                              ElementDofOrdering::LEXICOGRAPHIC,
-                              FaceType::Interior);
-
-   int n_int = R_int->Height();
-   Vector x(n_int), y1(n_int), y2(n_int);
-   x.Randomize(1);
-
-   // AddMultPA
-   y1 = 0.0;
-   y2 = 0.0;
-   integ1.AddMultPA(x, y1);
-   integ2.AddMultPA(x, y1);
-   integ_sum.AddMultPA(x, y2);
-   y1 -= y2;
-   REQUIRE(y1.Normlinf() == MFEM_Approx(0.0));
-
-   // AddMultTransposePA
-   y1 = 0.0;
-   y2 = 0.0;
-   integ1.AddMultTransposePA(x, y1);
-   integ2.AddMultTransposePA(x, y1);
-   integ_sum.AddMultTransposePA(x, y2);
-   y1 -= y2;
-   REQUIRE(y1.Normlinf() == MFEM_Approx(0.0));
-
-   // PA boundary
-   integ1.AssemblePABoundaryFaces(fes);
-   integ2.AssemblePABoundaryFaces(fes);
-   integ_sum.AssemblePABoundaryFaces(fes);
-
-   const Operator *R_bdr = fes.GetFaceRestriction(
-                              ElementDofOrdering::LEXICOGRAPHIC,
-                              FaceType::Boundary,
-                              L2FaceValues::DoubleValued);
-
-   int n_bdr = R_bdr->Height();
-   x.SetSize(n_bdr);
-   y1.SetSize(n_bdr);
-   y2.SetSize(n_bdr);
-   x.Randomize(1);
-
-   // AddMultPA
-   y1 = 0.0;
-   y2 = 0.0;
-   integ1.AddMultPA(x, y1);
-   integ2.AddMultPA(x, y1);
-   integ_sum.AddMultPA(x, y2);
-   y1 -= y2;
-   REQUIRE(y1.Normlinf() == MFEM_Approx(0.0));
-
-   // AddMultTransposePA
-   y1 = 0.0;
-   y2 = 0.0;
-   integ1.AddMultTransposePA(x, y1);
-   integ2.AddMultTransposePA(x, y1);
-   integ_sum.AddMultTransposePA(x, y2);
-   y1 -= y2;
-   REQUIRE(y1.Normlinf() == MFEM_Approx(0.0));
-}
-
-} // namespace pa_kernels