offset control

This yields better accuracy.

.gitignore

@@ -114,6 +114,9 @@ examples/ex28_*
 examples/ex28p_*
 examples/flux.*
 
+examples/euler-*
+examples/euler.mesh
+
 examples/amgx/ex1
 examples/amgx/ex1p
 examples/amgx/.logamgx

data/periodic-square-4x4.mesh

@@ -0,0 +1,127 @@
+MFEM mesh v1.0
+
+#
+# MFEM Geometry Types (see mesh/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+#
+
+dimension
+2
+
+# format: <attribute> <geometry type> <vertex 0> <vertex 1> ...
+elements
+16
+1  3 0 1 5 4
+2  3 1 2 6 5
+3  3 2 3 7 6
+4  3 3 0 4 7
+5  3 4 5 9 8
+6  3 5 6 10 9
+7  3 6 7 11 10
+8  3 7 4 8 11
+9  3 8 9 13 12
+10 3 9 10 14 13
+11 3 10 11 15 14
+12 3 11 8 12 15
+13 3 12 13  1  0
+14 3 13 14  2  1
+15 3 14 15  3  2
+16 3 15 12  0  3
+
+boundary
+0
+
+vertices
+16
+
+nodes
+FiniteElementSpace
+FiniteElementCollection: L2_T1_2D_P1
+VDim: 2
+Ordering: 1
+
+-1.0 -1.0
+-0.5 -1.0
+-1.0 -0.5
+-0.5 -0.5
+
+-0.5 -1.0
++0.0 -1.0
+-0.5 -0.5
++0.0 -0.5
+
++0.0 -1.0
++0.5 -1.0
++0.0 -0.5
++0.5 -0.5
+
++0.5 -1.0
++1.0 -1.0
++0.5 -0.5
++1.0 -0.5
+
+-1.0 -0.5
+-0.5 -0.5
+-1.0 +0.0
+-0.5 +0.0
+
+-0.5 -0.5
++0.0 -0.5
+-0.5 +0.0
++0.0 +0.0
+
++0.0 -0.5
++0.5 -0.5
++0.0 +0.0
++0.5 +0.0
+
++0.5 -0.5
++1.0 -0.5
++0.5 +0.0
++1.0 +0.0
+
+-1.0 +0.0
+-0.5 +0.0
+-1.0 +0.5
+-0.5 +0.5
+
+-0.5 +0.0
++0.0 +0.0
+-0.5 +0.5
++0.0 +0.5
+
++0.0 +0.0
++0.5 +0.0
++0.0 +0.5
++0.5 +0.5
+
++0.5 +0.0
++1.0 +0.0
++0.5 +0.5
++1.0 +0.5
+
+-1.0 +0.5
+-0.5 +0.5
+-1.0 +1.0
+-0.5 +1.0
+
+-0.5 +0.5
++0.0 +0.5
+-0.5 +1.0
++0.0 +1.0
+
++0.0 +0.5
++0.5 +0.5
++0.0 +1.0
++0.5 +1.0
+
++0.5 +0.5
++1.0 +0.5
++0.5 +1.0
++1.0 +1.0

examples/adv.cpp

@@ -0,0 +1,303 @@
+// Linear advection uₜ + ∇⋅F(u) = 0
+// with bound constraint, u ≥ uₘ
+//
+//
+
+#include <fstream>
+#include <iostream>
+#include <sstream>
+#include <cmath>
+
+#include "mfem.hpp"
+
+// Classes HyperbolicConservationLaws, RiemannSolver, and FaceIntegrator
+// shared between the serial and parallel version of the example.
+#include "proxGalerkinHCL.hpp"
+
+using namespace std;
+using namespace mfem;
+
+void AdvMesh(const int problem, const char **mesh_file);
+
+VectorFunctionCoefficient AdvInitCondition(const int problem);
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   int problem = 1;
+
+   const char *mesh_file = "";
+   int intOrderOffset = 3;
+   int ref_levels = 2;
+   int order = 3;
+   int ode_solver_type = 11;
+   double t_final = 2.0;
+   double dt = -0.01;
+   double cfl = 0.3;
+   bool visualization = true;
+   int vis_steps = 50;
+
+   int precision = 8;
+   out.precision(precision);
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
+   args.AddOption(&problem, "-p", "--problem",
+                  "Problem setup to use. See options in velocity_function().");
+   args.AddOption(&ref_levels, "-r", "--refine",
+                  "Number of times to refine the mesh uniformly.");
+   args.AddOption(&order, "-o", "--order",
+                  "Order (degree) of the finite elements.");
+   args.AddOption(&ode_solver_type, "-s", "--ode-solver",
+                  "ODE solver: 1 - Forward Euler,\n\t"
+                  "            2 - RK2 SSP, 3 - RK3 SSP, 4 - RK4, 6 - RK6.");
+   args.AddOption(&t_final, "-tf", "--t-final", "Final time; start time is 0.");
+   args.AddOption(&dt, "-dt", "--time-step",
+                  "Time step. Positive number skips CFL timestep calculation.");
+   args.AddOption(&cfl, "-c", "--cfl-number",
+                  "CFL number for timestep calculation.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&vis_steps, "-vs", "--visualization-steps",
+                  "Visualize every n-th timestep.");
+
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(out);
+      return 1;
+   }
+   // When the user does not provide mesh file,
+   // use the default mesh file for the problem.
+   if ((mesh_file == NULL) || (mesh_file[0] == '\0'))    // if NULL or empty
+   {
+      AdvMesh(problem, &mesh_file);  // get default mesh file name
+   }
+
+   // 2. Read the mesh from the given mesh file.
+   Mesh mesh = Mesh(mesh_file);
+   const int dim = mesh.Dimension();
+   const int num_equations = 1;
+
+   if (problem == 5)
+   {
+      mesh.Transform([](const Vector &x, Vector &y)
+      {
+         y = x;
+         y *= 0.5;
+      });
+   }
+   // perform uniform refine
+   for (int lev = 0; lev < ref_levels; lev++)
+   {
+      mesh.UniformRefinement();
+   }
+   if (dim > 1) { mesh.EnsureNCMesh(); }
+
+   // 3. Define the ODE solver used for time integration. Several explicit
+   //    Runge-Kutta methods are available.
+   ODESolver *ode_solver = NULL;
+   switch (ode_solver_type)
+   {
+      case 11:
+         ode_solver = new BackwardEulerSolver;
+         break;
+      default:
+         out << "Unknown ODE solver type: " << ode_solver_type << '\n';
+         return 3;
+   }
+
+   // 4. Define the discontinuous DG finite element space of the given
+   //    polynomial order on the refined mesh.
+   DG_FECollection fec(order, dim);
+   // Finite element space for a scalar (thermodynamic quantity)
+   FiniteElementSpace fes(&mesh, &fec);
+   // Finite element space for a mesh-dim vector quantity (momentum)
+   FiniteElementSpace dfes(&mesh, &fec, dim, Ordering::byNODES);
+   // Finite element space for all variables together (total thermodynamic state)
+   FiniteElementSpace vfes(&mesh, &fec, num_equations, Ordering::byNODES);
+
+   // This example depends on this ordering of the space.
+   out << "Number of unknowns: " << vfes.GetVSize() << endl;
+
+   // 6. Define the initial conditions, save the corresponding mesh and grid
+   //    functions to a file. This can be opened with GLVis with the -gc option.
+   // Initialize the state.
+   VectorFunctionCoefficient u0 = AdvInitCondition(problem);
+   GridFunction sol(&vfes);
+   sol.ProjectCoefficient(u0);
+
+   // Output the initial solution.
+   {
+      ostringstream mesh_name;
+      mesh_name << "adv-mesh.mesh";
+      ofstream mesh_ofs(mesh_name.str().c_str());
+      mesh_ofs.precision(precision);
+      mesh_ofs << mesh;
+
+      for (int k = 0; k < num_equations; k++)
+      {
+         GridFunction uk(&fes, sol.GetData() + k * fes.GetNDofs());
+         ostringstream sol_name;
+         sol_name << "adv-" << k << "-init.gf";
+         ofstream sol_ofs(sol_name.str().c_str());
+         sol_ofs.precision(precision);
+         sol_ofs << uk;
+      }
+   }
+
+   // 7. Set up the nonlinear form corresponding to the DG discretization of the
+   //    flux divergence, and assemble the corresponding mass matrix.
+   RiemannSolver *numericalFlux = new RusanovFlux();
+   Vector b(2); b[0] = 1.0; b[1] = 1.0;
+   VectorConstantCoefficient b_cf(b);
+   AdvectionFormIntegrator adv_fi(new RusanovFlux(), dim, b_cf, intOrderOffset);
+   ProxGalerkinPolynomialAlphaMaker alphamaker(1.0, 1.0);
+   ProxGalerkinHCL adv(&fes, &vfes, adv_fi, num_equations, alphamaker);
+
+   // Visualize the density
+   socketstream sout;
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int visport = 19916;
+
+      sout.open(vishost, visport);
+      if (!sout)
+      {
+         visualization = false;
+         out << "Unable to connect to GLVis server at " << vishost << ':'
+             << visport << endl;
+         out << "GLVis visualization disabled.\n";
+      }
+      else
+      {
+         GridFunction mom(&dfes, sol.GetData());
+         sout << "solution\n" << mesh << mom;
+         sout << "view 0 0\n";  // view from top
+         sout << "keys jlm\n";  // turn off perspective and light
+         sout << "pause\n";
+         sout << flush;
+         out << "GLVis visualization paused."
+             << " Press space (in the GLVis window) to resume it.\n";
+      }
+   }
+
+   // Determine the minimum element size.
+   double hmin = infinity();
+   if (cfl > 0)
+   {
+      for (int i = 0; i < mesh.GetNE(); i++)
+      {
+         hmin = min(mesh.GetElementSize(i, 1), hmin);
+      }
+   }
+
+   // Start the timer.
+   tic_toc.Clear();
+   tic_toc.Start();
+
+   double t = 0.0;
+   adv.SetTime(t);
+   ode_solver->Init(adv);
+
+   if (cfl > 0)
+   {
+      // Find a safe dt, using a temporary vector. Calling Mult() computes the
+      // maximum char speed at all quadrature points on all faces.
+      Vector z(sol.Size());
+      adv.Mult(sol, z);
+
+      double max_char_speed = adv.getMaxCharSpeed();
+      dt = cfl * hmin / max_char_speed / (2 * order + 1);
+   }
+
+   // Integrate in time.
+   bool done = false;
+   for (int ti = 0; !done;)
+   {
+      double dt_real = min(dt, t_final - t);
+
+      ode_solver->Step(sol, t, dt_real);
+      if (cfl > 0)
+      {
+         double max_char_speed = adv.getMaxCharSpeed();
+         dt = cfl * hmin / max_char_speed / (2 * order + 1);
+      }
+      ti++;
+
+      done = (t >= t_final - 1e-8 * dt);
+      if (done || ti % vis_steps == 0)
+      {
+         out << "time step: " << ti << ", time: " << t << endl;
+         if (visualization)
+         {
+            GridFunction mom(&dfes, sol.GetData());
+            sout << "solution\n" << mesh << mom << flush;
+            sout << "window_title 't = " << t << "'";
+         }
+      }
+   }
+
+   tic_toc.Stop();
+   out << " done, " << tic_toc.RealTime() << "s." << endl;
+
+   // 9. Save the final solution. This output can be viewed later using GLVis:
+   //    "glvis -m adv.mesh -g adv-1-final.gf".
+   {
+      ostringstream mesh_name;
+      mesh_name << "adv-mesh-final.mesh";
+      ofstream mesh_ofs(mesh_name.str().c_str());
+      mesh_ofs.precision(precision);
+      mesh_ofs << mesh;
+
+      for (int k = 0; k < num_equations; k++)
+      {
+         GridFunction uk(&fes, sol.GetData() + k * fes.GetNDofs());
+         ostringstream sol_name;
+         sol_name << "adv-" << k << "-final.gf";
+         ofstream sol_ofs(sol_name.str().c_str());
+         sol_ofs.precision(precision);
+         sol_ofs << uk;
+      }
+   }
+
+   // 10. Compute the L2 solution error summed for all components.
+   //   if (t_final == 2.0) {
+   const double error = sol.ComputeLpError(2, u0);
+   out << "Solution error: " << error << endl;
+
+   // Free the used memory.
+   delete ode_solver;
+
+   return 0;
+}
+
+void AdvMesh(const int problem, const char **mesh_file)
+{
+   switch (problem)
+   {
+      case 1:
+         *mesh_file = "../data/periodic-square-4x4.mesh";
+         break;
+      default:
+         throw invalid_argument("Default mesh is undefined");
+   }
+}
+
+// Initial condition
+VectorFunctionCoefficient AdvInitCondition(const int problem)
+{
+   switch (problem)
+   {
+      case 1: // fast moving vortex
+         return VectorFunctionCoefficient(1, [](const Vector &x, Vector &y)
+         {
+            MFEM_ASSERT(x.Size() == 2, "");
+            y = std::sin(2.0*M_PI*x[0])*std::sin(2.0*M_PI*x[1]) + 1.0;
+         });
+      default:
+         throw invalid_argument("Problem Undefined");
+   }
+}

examples/draw_solutions.m

@@ -0,0 +1,31 @@
+meshfiles = dir('mesh*.mesh');
+ufiles = dir('u*.gf');
+rhofiles = dir('rho*.gf');
+frhofiles = dir('f_rho*.gf');
+
+ax = arrayfun(@(i) axes(figure()), 1:3);
+figs = arrayfun(@(ax) ax.Parent, ax);
+recorders = arrayfun(@(name, fig) FigureRecorder(name, fig), ["u.mp4", "rho.mp4", "frho.mp4"], figs);
+for i = 1 : length(meshfiles)
+    mesh = read_mfem_mesh([meshfiles(i).folder filesep meshfiles(i).name]);
+
+    contour_mfem_gf(ax(1), mesh, [ufiles(i).folder filesep ufiles(i).name], 11);
+    title(ax(1), sprintf('Iteration %d: u', i))
+
+    contour_mfem_gf(ax(2), mesh, [rhofiles(i).folder filesep rhofiles(i).name], 0.1:0.1:0.9);
+    title(ax(2), sprintf('Iteration %d: ρ', i))
+    set(ax(2), 'Clim', [0, 1]);
+    colormap(ax(2), flipud(gray));
+
+    contour_mfem_gf(ax(3), mesh, [frhofiles(i).folder filesep frhofiles(i).name], 0.1:0.1:0.9);
+    title(ax(3), sprintf('Iteration %d: ρ̃', i))
+    set(ax(3), 'Clim', [0, 1]);
+    colormap(ax(3), flipud(gray));
+
+    arrayfun(@(a) view(a, 2), ax);
+    arrayfun(@(a) set(a, 'FontSize', 16), ax);
+    arrayfun(@(a) colorbar(a), ax);
+    arrayfun(@(a) axis(a, 'off'), ax)
+    arrayfun(@(recorder) recorder.capture(), recorders);
+end
+arrayfun(@(recorder) recorder.close(), recorders);

examples/ex18.cpp

@@ -1,6 +1,6 @@
 //                                MFEM Example 18
 //
-// Compile with: make ex18
+// Compile with: make ex18p
 //
 // Sample runs:
 //
@@ -14,6 +14,8 @@
 //               equations, a model nonlinear hyperbolic PDE, with a
 //               discontinuous Galerkin (DG) formulation.
 //
+//                (u_t, v)_T - (F(u), ∇ v)_T + (F̂(u,n), [[v]])_F = 0
+//
 //               Specifically, it solves for an exact solution of the equations
 //               whereby a vortex is transported by a uniform flow. Since all
 //               boundaries are periodic here, the method's accuracy can be
@@ -27,43 +29,46 @@
 //               method. An additional factor can be tuned by passing the --cfl
 //               (or -c shorter) flag.
 //
-//               The example demonstrates user-defined bilinear and nonlinear
-//               form integrators for systems of equations that are defined with
-//               block vectors, and how these are used with an operator for
-//               explicit time integrators. In this case the system also
-//               involves an external approximate Riemann solver for the DG
-//               interface flux. It also demonstrates how to use GLVis for
-//               in-situ visualization of vector grid functions.
+//               The example demonstrates usage of DGHyperbolicConservationLaws
+//               that wraps NonlinearFormIntegrators containing element and face
+//               integration schemes. In this case the system also involves an
+//               external approximate Riemann solver for the DG interface flux.
+//               It also demonstrates how to use GLVis for in-situ visualization
+//               of vector grid function and how to set top-view.
 //
 //               We recommend viewing examples 9, 14 and 17 before viewing this
 //               example.
 
-#include "mfem.hpp"
 #include <fstream>
-#include <sstream>
 #include <iostream>
+#include <sstream>
+#include <cmath>
 
-// Classes FE_Evolution, RiemannSolver, and FaceIntegrator
+#include "mfem.hpp"
+
+// Classes HyperbolicConservationLaws, RiemannSolver, and FaceIntegrator
 // shared between the serial and parallel version of the example.
-#include "ex18.hpp"
+#include "fem/hyperbolic_conservation_laws.hpp"
 
-// Choice for the problem setup. See InitialCondition in ex18.hpp.
-int problem;
+using namespace std;
+using namespace mfem;
 
-// Equation constant parameters.
-const int num_equation = 4;
-const double specific_heat_ratio = 1.4;
-const double gas_constant = 1.0;
+void EulerMesh(const int problem, const char **mesh_file);
 
-// Maximum characteristic speed (updated by integrators)
-double max_char_speed;
+VectorFunctionCoefficient EulerInitialCondition(const int problem,
+                                                const double specific_heat_ratio,
+                                                const double gas_constant);
 
 int main(int argc, char *argv[])
 {
    // 1. Parse command-line options.
-   problem = 1;
-   const char *mesh_file = "../data/periodic-square.mesh";
-   int ref_levels = 1;
+   int problem = 1;
+   const double specific_heat_ratio = 1.4;
+   const double gas_constant = 1.0;
+
+   const char *mesh_file = "";
+   int IntOrderOffset = 3;
+   int ref_levels = 2;
    int order = 3;
    int ode_solver_type = 4;
    double t_final = 2.0;
@@ -73,11 +78,10 @@ int main(int argc, char *argv[])
    int vis_steps = 50;
 
    int precision = 8;
-   cout.precision(precision);
+   out.precision(precision);
 
    OptionsParser args(argc, argv);
-   args.AddOption(&mesh_file, "-m", "--mesh",
-                  "Mesh file to use.");
+   args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
    args.AddOption(&problem, "-p", "--problem",
                   "Problem setup to use. See options in velocity_function().");
    args.AddOption(&ref_levels, "-r", "--refine",
@@ -87,8 +91,7 @@ int main(int argc, char *argv[])
    args.AddOption(&ode_solver_type, "-s", "--ode-solver",
                   "ODE solver: 1 - Forward Euler,\n\t"
                   "            2 - RK2 SSP, 3 - RK3 SSP, 4 - RK4, 6 - RK6.");
-   args.AddOption(&t_final, "-tf", "--t-final",
-                  "Final time; start time is 0.");
+   args.AddOption(&t_final, "-tf", "--t-final", "Final time; start time is 0.");
    args.AddOption(&dt, "-dt", "--time-step",
                   "Time step. Positive number skips CFL timestep calculation.");
    args.AddOption(&cfl, "-c", "--cfl-number",
@@ -102,42 +105,62 @@ int main(int argc, char *argv[])
    args.Parse();
    if (!args.Good())
    {
-      args.PrintUsage(cout);
+      args.PrintUsage(out);
       return 1;
    }
-   args.PrintOptions(cout);
+   // When the user does not provide mesh file,
+   // use the default mesh file for the problem.
+   if ((mesh_file == NULL) || (mesh_file[0] == '\0'))    // if NULL or empty
+   {
+      EulerMesh(problem, &mesh_file);  // get default mesh file name
+   }
 
-   // 2. Read the mesh from the given mesh file. This example requires a 2D
-   //    periodic mesh, such as ../data/periodic-square.mesh.
-   Mesh mesh(mesh_file, 1, 1);
+   // 2. Read the mesh from the given mesh file.
+   Mesh mesh = Mesh(mesh_file);
    const int dim = mesh.Dimension();
+   const int num_equations = dim + 2;
 
-   MFEM_ASSERT(dim == 2, "Need a two-dimensional mesh for the problem definition");
+   if (problem == 5)
+   {
+      mesh.Transform([](const Vector &x, Vector &y)
+      {
+         y = x;
+         y *= 0.5;
+      });
+   }
+   // perform uniform refine
+   for (int lev = 0; lev < ref_levels; lev++)
+   {
+      mesh.UniformRefinement();
+   }
+   if (dim > 1) { mesh.EnsureNCMesh(); }
 
    // 3. Define the ODE solver used for time integration. Several explicit
    //    Runge-Kutta methods are available.
    ODESolver *ode_solver = NULL;
    switch (ode_solver_type)
    {
-      case 1: ode_solver = new ForwardEulerSolver; break;
-      case 2: ode_solver = new RK2Solver(1.0); break;
-      case 3: ode_solver = new RK3SSPSolver; break;
-      case 4: ode_solver = new RK4Solver; break;
-      case 6: ode_solver = new RK6Solver; break;
+      case 1:
+         ode_solver = new ForwardEulerSolver;
+         break;
+      case 2:
+         ode_solver = new RK2Solver(1.0);
+         break;
+      case 3:
+         ode_solver = new RK3SSPSolver;
+         break;
+      case 4:
+         ode_solver = new RK4Solver;
+         break;
+      case 6:
+         ode_solver = new RK6Solver;
+         break;
       default:
-         cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
+         out << "Unknown ODE solver type: " << ode_solver_type << '\n';
          return 3;
    }
 
-   // 4. Refine the mesh to increase the resolution. In this example we do
-   //    'ref_levels' of uniform refinement, where 'ref_levels' is a
-   //    command-line parameter.
-   for (int lev = 0; lev < ref_levels; lev++)
-   {
-      mesh.UniformRefinement();
-   }
-
-   // 5. Define the discontinuous DG finite element space of the given
+   // 4. Define the discontinuous DG finite element space of the given
    //    polynomial order on the refined mesh.
    DG_FECollection fec(order, dim);
    // Finite element space for a scalar (thermodynamic quantity)
@@ -145,41 +168,34 @@ int main(int argc, char *argv[])
    // Finite element space for a mesh-dim vector quantity (momentum)
    FiniteElementSpace dfes(&mesh, &fec, dim, Ordering::byNODES);
    // Finite element space for all variables together (total thermodynamic state)
-   FiniteElementSpace vfes(&mesh, &fec, num_equation, Ordering::byNODES);
+   FiniteElementSpace vfes(&mesh, &fec, num_equations, Ordering::byNODES);
 
    // This example depends on this ordering of the space.
    MFEM_ASSERT(fes.GetOrdering() == Ordering::byNODES, "");
 
-   cout << "Number of unknowns: " << vfes.GetVSize() << endl;
+   out << "Number of unknowns: " << vfes.GetVSize() << endl;
 
    // 6. Define the initial conditions, save the corresponding mesh and grid
    //    functions to a file. This can be opened with GLVis with the -gc option.
-
-   // The solution u has components {density, x-momentum, y-momentum, energy}.
-   // These are stored contiguously in the BlockVector u_block.
-   Array<int> offsets(num_equation + 1);
-   for (int k = 0; k <= num_equation; k++) { offsets[k] = k * vfes.GetNDofs(); }
-   BlockVector u_block(offsets);
-
-   // Momentum grid function on dfes for visualization.
-   GridFunction mom(&dfes, u_block.GetData() + offsets[1]);
-
    // Initialize the state.
-   VectorFunctionCoefficient u0(num_equation, InitialCondition);
-   GridFunction sol(&vfes, u_block.GetData());
+   VectorFunctionCoefficient u0 = EulerInitialCondition(problem,
+                                                        specific_heat_ratio, gas_constant);
+   GridFunction sol(&vfes);
    sol.ProjectCoefficient(u0);
 
    // Output the initial solution.
    {
-      ofstream mesh_ofs("vortex.mesh");
+      ostringstream mesh_name;
+      mesh_name << "euler-mesh.mesh";
+      ofstream mesh_ofs(mesh_name.str().c_str());
       mesh_ofs.precision(precision);
       mesh_ofs << mesh;
 
-      for (int k = 0; k < num_equation; k++)
+      for (int k = 0; k < num_equations; k++)
       {
-         GridFunction uk(&fes, u_block.GetBlock(k));
+         GridFunction uk(&fes, sol.GetData() + k * fes.GetNDofs());
          ostringstream sol_name;
-         sol_name << "vortex-" << k << "-init.gf";
+         sol_name << "euler-" << k << "-init.gf";
          ofstream sol_ofs(sol_name.str().c_str());
          sol_ofs.precision(precision);
          sol_ofs << uk;
@@ -188,51 +204,43 @@ int main(int argc, char *argv[])
 
    // 7. Set up the nonlinear form corresponding to the DG discretization of the
    //    flux divergence, and assemble the corresponding mass matrix.
-   MixedBilinearForm Aflux(&dfes, &fes);
-   Aflux.AddDomainIntegrator(new TransposeIntegrator(new GradientIntegrator()));
-   Aflux.Assemble();
-
-   NonlinearForm A(&vfes);
-   RiemannSolver rsolver;
-   A.AddInteriorFaceIntegrator(new FaceIntegrator(rsolver, dim));
-
-   // 8. 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 euler(vfes, A, Aflux.SpMat());
+   RiemannSolver *numericalFlux = new RusanovFlux();
+   DGHyperbolicConservationLaws euler = getEulerSystem(
+                                           &vfes, numericalFlux, specific_heat_ratio, IntOrderOffset);
 
    // Visualize the density
    socketstream sout;
    if (visualization)
    {
       char vishost[] = "localhost";
-      int  visport   = 19916;
+      int visport = 19916;
 
       sout.open(vishost, visport);
       if (!sout)
       {
-         cout << "Unable to connect to GLVis server at "
-              << vishost << ':' << visport << endl;
          visualization = false;
-         cout << "GLVis visualization disabled.\n";
+         out << "Unable to connect to GLVis server at " << vishost << ':'
+             << visport << endl;
+         out << "GLVis visualization disabled.\n";
       }
       else
       {
-         sout.precision(precision);
+         GridFunction mom(&dfes, sol.GetData());
          sout << "solution\n" << mesh << mom;
+         sout << "view 0 0\n";  // view from top
+         sout << "keys jlm\n";  // turn off perspective and light
          sout << "pause\n";
          sout << flush;
-         cout << "GLVis visualization paused."
-              << " Press space (in the GLVis window) to resume it.\n";
+         out << "GLVis visualization paused."
+             << " Press space (in the GLVis window) to resume it.\n";
       }
    }
 
    // Determine the minimum element size.
-   double hmin = 0.0;
+   double hmin = infinity();
    if (cfl > 0)
    {
-      hmin = mesh.GetElementSize(0, 1);
-      for (int i = 1; i < mesh.GetNE(); i++)
+      for (int i = 0; i < mesh.GetNE(); i++)
       {
          hmin = min(mesh.GetElementSize(i, 1), hmin);
       }
@@ -250,60 +258,258 @@ int main(int argc, char *argv[])
    {
       // Find a safe dt, using a temporary vector. Calling Mult() computes the
       // maximum char speed at all quadrature points on all faces.
-      Vector z(A.Width());
-      max_char_speed = 0.;
-      A.Mult(sol, z);
-      dt = cfl * hmin / max_char_speed / (2*order+1);
+      Vector z(sol.Size());
+      euler.Mult(sol, z);
+
+      double max_char_speed = euler.getMaxCharSpeed();
+      dt = cfl * hmin / max_char_speed / (2 * order + 1);
    }
 
    // Integrate in time.
    bool done = false;
-   for (int ti = 0; !done; )
+   for (int ti = 0; !done;)
    {
       double dt_real = min(dt, t_final - t);
 
       ode_solver->Step(sol, t, dt_real);
       if (cfl > 0)
       {
-         dt = cfl * hmin / max_char_speed / (2*order+1);
+         double max_char_speed = euler.getMaxCharSpeed();
+         dt = cfl * hmin / max_char_speed / (2 * order + 1);
       }
       ti++;
 
-      done = (t >= t_final - 1e-8*dt);
+      done = (t >= t_final - 1e-8 * dt);
       if (done || ti % vis_steps == 0)
       {
-         cout << "time step: " << ti << ", time: " << t << endl;
+         out << "time step: " << ti << ", time: " << t << endl;
          if (visualization)
          {
+            GridFunction mom(&dfes, sol.GetData());
             sout << "solution\n" << mesh << mom << flush;
+            sout << "window_title 't = " << t << "'";
          }
       }
    }
 
    tic_toc.Stop();
-   cout << " done, " << tic_toc.RealTime() << "s." << endl;
+   out << " done, " << tic_toc.RealTime() << "s." << endl;
 
    // 9. Save the final solution. This output can be viewed later using GLVis:
-   //    "glvis -m vortex.mesh -g vortex-1-final.gf".
-   for (int k = 0; k < num_equation; k++)
+   //    "glvis -m euler.mesh -g euler-1-final.gf".
    {
-      GridFunction uk(&fes, u_block.GetBlock(k));
-      ostringstream sol_name;
-      sol_name << "vortex-" << k << "-final.gf";
-      ofstream sol_ofs(sol_name.str().c_str());
-      sol_ofs.precision(precision);
-      sol_ofs << uk;
+      ostringstream mesh_name;
+      mesh_name << "euler-mesh-final.mesh";
+      ofstream mesh_ofs(mesh_name.str().c_str());
+      mesh_ofs.precision(precision);
+      mesh_ofs << mesh;
+
+      for (int k = 0; k < num_equations; k++)
+      {
+         GridFunction uk(&fes, sol.GetData() + k * fes.GetNDofs());
+         ostringstream sol_name;
+         sol_name << "euler-" << k << "-final.gf";
+         ofstream sol_ofs(sol_name.str().c_str());
+         sol_ofs.precision(precision);
+         sol_ofs << uk;
+      }
    }
 
    // 10. Compute the L2 solution error summed for all components.
-   if (t_final == 2.0)
-   {
-      const double error = sol.ComputeLpError(2, u0);
-      cout << "Solution error: " << error << endl;
-   }
+   //   if (t_final == 2.0) {
+   const double error = sol.ComputeLpError(2, u0);
+   out << "Solution error: " << error << endl;
 
    // Free the used memory.
    delete ode_solver;
 
    return 0;
 }
+
+void EulerMesh(const int problem, const char **mesh_file)
+{
+   switch (problem)
+   {
+      case 1:
+         *mesh_file = "../data/periodic-square-4x4.mesh";
+         break;
+      case 2:
+         *mesh_file = "../data/periodic-square-4x4.mesh";
+         break;
+      case 3:
+         *mesh_file = "../data/periodic-square-4x4.mesh";
+         break;
+      case 4:
+         *mesh_file = "../data/periodic-segment.mesh";
+         break;
+      case 5:
+         *mesh_file = "../data/periodic-square-4x4.mesh";
+         break;
+      default:
+         throw invalid_argument("Default mesh is undefined");
+   }
+}
+
+// Initial condition
+VectorFunctionCoefficient EulerInitialCondition(const int problem,
+                                                const double specific_heat_ratio,
+                                                const double gas_constant)
+{
+   switch (problem)
+   {
+      case 1: // fast moving vortex
+         return VectorFunctionCoefficient(4, [specific_heat_ratio,
+                                              gas_constant](const Vector &x, Vector &y)
+         {
+            MFEM_ASSERT(x.Size() == 2, "");
+
+            double radius = 0, Minf = 0, beta = 0;
+            // "Fast euler"
+            radius = 0.2;
+            Minf = 0.5;
+            beta = 1. / 5.;
+
+            const double xc = 0.0, yc = 0.0;
+
+            // Nice units
+            const double vel_inf = 1.;
+            const double den_inf = 1.;
+
+            // Derive remainder of background state from this and Minf
+            const double pres_inf = (den_inf / specific_heat_ratio) *
+                                    (vel_inf / Minf) * (vel_inf / Minf);
+            const double temp_inf = pres_inf / (den_inf * gas_constant);
+
+            double r2rad = 0.0;
+            r2rad += (x(0) - xc) * (x(0) - xc);
+            r2rad += (x(1) - yc) * (x(1) - yc);
+            r2rad /= (radius * radius);
+
+            const double shrinv1 = 1.0 / (specific_heat_ratio - 1.);
+
+            const double velX =
+               vel_inf * (1 - beta * (x(1) - yc) / radius * exp(-0.5 * r2rad));
+            const double velY =
+               vel_inf * beta * (x(0) - xc) / radius * exp(-0.5 * r2rad);
+            const double vel2 = velX * velX + velY * velY;
+
+            const double specific_heat =
+               gas_constant * specific_heat_ratio * shrinv1;
+            const double temp = temp_inf - 0.5 * (vel_inf * beta) *
+                                (vel_inf * beta) / specific_heat *
+                                exp(-r2rad);
+
+            const double den = den_inf * pow(temp / temp_inf, shrinv1);
+            const double pres = den * gas_constant * temp;
+            const double energy = shrinv1 * pres / den + 0.5 * vel2;
+
+            y(0) = den;
+            y(1) = den * velX;
+            y(2) = den * velY;
+            y(3) = den * energy;
+         });
+      case 2: // slow moving vortex
+         return VectorFunctionCoefficient(4, [specific_heat_ratio,
+                                              gas_constant](const Vector &x, Vector &y)
+         {
+            MFEM_ASSERT(x.Size() == 2, "");
+
+            double radius = 0, Minf = 0, beta = 0;
+            // "Slow euler"
+            radius = 0.2;
+            Minf = 0.05;
+            beta = 1. / 50.;
+
+            const double xc = 0.0, yc = 0.0;
+
+            // Nice units
+            const double vel_inf = 1.;
+            const double den_inf = 1.;
+
+            // Derive remainder of background state from this and Minf
+            const double pres_inf = (den_inf / specific_heat_ratio) *
+                                    (vel_inf / Minf) * (vel_inf / Minf);
+            const double temp_inf = pres_inf / (den_inf * gas_constant);
+
+            double r2rad = 0.0;
+            r2rad += (x(0) - xc) * (x(0) - xc);
+            r2rad += (x(1) - yc) * (x(1) - yc);
+            r2rad /= (radius * radius);
+
+            const double shrinv1 = 1.0 / (specific_heat_ratio - 1.);
+
+            const double velX =
+               vel_inf * (1 - beta * (x(1) - yc) / radius * exp(-0.5 * r2rad));
+            const double velY =
+               vel_inf * beta * (x(0) - xc) / radius * exp(-0.5 * r2rad);
+            const double vel2 = velX * velX + velY * velY;
+
+            const double specific_heat =
+               gas_constant * specific_heat_ratio * shrinv1;
+            const double temp = temp_inf - 0.5 * (vel_inf * beta) *
+                                (vel_inf * beta) / specific_heat *
+                                exp(-r2rad);
+
+            const double den = den_inf * pow(temp / temp_inf, shrinv1);
+            const double pres = den * gas_constant * temp;
+            const double energy = shrinv1 * pres / den + 0.5 * vel2;
+
+            y(0) = den;
+            y(1) = den * velX;
+            y(2) = den * velY;
+            y(3) = den * energy;
+         });
+      case 3: // moving sine wave
+         return VectorFunctionCoefficient(4, [](const Vector &x, Vector &y)
+         {
+            MFEM_ASSERT(x.Size() == 2, "");
+            const double density = 1.0 + 0.2 * sin(M_PI*(x(0) + x(1)));
+            const double velocity_x = 0.7;
+            const double velocity_y = 0.3;
+            const double pressure = 1.0;
+            const double energy =
+               pressure / (1.4 - 1.0) +
+               density * 0.5 * (velocity_x * velocity_x + velocity_y * velocity_y);
+
+            y(0) = density;
+            y(1) = density * velocity_x;
+            y(2) = density * velocity_y;
+            y(3) = energy;
+         });
+      case 4:
+         return VectorFunctionCoefficient(3, [](const Vector &x, Vector &y)
+         {
+            MFEM_ASSERT(x.Size() == 1, "");
+            const double density = 1.0 + 0.2 * sin(M_PI * 2 * x(0));
+            const double velocity_x = 1.0;
+            const double pressure = 1.0;
+            const double energy =
+               pressure / (1.4 - 1.0) + density * 0.5 * (velocity_x * velocity_x);
+
+            y(0) = density;
+            y(1) = density * velocity_x;
+            y(2) = energy;
+         });
+      case 5:
+         return VectorFunctionCoefficient(4, [](const Vector &x, Vector &y)
+         {
+            MFEM_ASSERT(x.Size() == 2, "");
+            const double L = 1.0;
+            const double density = abs(x(1)) < 0.25 ? 2 : 1;
+            const double velocity_x = abs(x(1)) < 0.25 ? -0.5 : 0.5;
+            const double velocity_y = abs(x(1)) < 0.25 ? 0.01 * sin(M_PI*x(0) / L)
+                                      : 0.01 * sin(M_PI*x(0) / L);
+            const double pressure = abs(x(1)) < 0.25 ? 2.5 : 2.5;
+            const double energy =
+               pressure / (1.4 - 1.0) +
+               density * 0.5 * (velocity_x * velocity_x + velocity_y * velocity_y);
+
+            y(0) = density;
+            y(1) = density * velocity_x;
+            y(2) = density * velocity_y;
+            y(3) = energy;
+         });
+      default:
+         throw invalid_argument("Problem Undefined");
+   }
+}

examples/ex18p.cpp

@@ -1,24 +1,26 @@
-//                       MFEM Example 18 - Parallel Version
+//                                MFEM Example 18 - Parallel Version
 //
-// Compile with: make ex18
+// Compile with: make ex18p
 //
 // Sample runs:
 //
-//       mpirun -np 4 ex18p -p 1 -rs 2 -rp 1 -o 1 -s 3
-//       mpirun -np 4 ex18p -p 1 -rs 1 -rp 1 -o 3 -s 4
-//       mpirun -np 4 ex18p -p 1 -rs 1 -rp 1 -o 5 -s 6
-//       mpirun -np 4 ex18p -p 2 -rs 1 -rp 1 -o 1 -s 3
-//       mpirun -np 4 ex18p -p 2 -rs 1 -rp 1 -o 3 -s 3
+//       mpirun ex18p -p 1 -rs 2 -o 1 -s 3
+//       mpirun ex18p -p 1 -rs 1 -o 3 -s 4
+//       mpirun ex18p -p 1 -rs 0 -o 5 -s 6
+//       mpirun ex18p -p 2 -rs 1 -o 1 -s 3
+//       mpirun ex18p -p 2 -rs 0 -o 3 -s 3
 //
 // Description:  This example code solves the compressible Euler system of
 //               equations, a model nonlinear hyperbolic PDE, with a
-//               discontinuous Galerkin (DG) formulation.
+//               discontinuous Galerkin (DG) formulation in parallel.
+//
+//                (u_t, v)_T - (F(u), ∇ v)_T + (F̂(u,n), [[v]])_F = 0
 //
 //               Specifically, it solves for an exact solution of the equations
-//               whereby a vortex is transported by a uniform flow. Since all
+//               whereby a euler is transported by a uniform flow. Since all
 //               boundaries are periodic here, the method's accuracy can be
 //               assessed by measuring the difference between the solution and
-//               the initial condition at a later time when the vortex returns
+//               the initial condition at a later time when the euler returns
 //               to its initial location.
 //
 //               Note that as the order of the spatial discretization increases,
@@ -27,48 +29,52 @@
 //               method. An additional factor can be tuned by passing the --cfl
 //               (or -c shorter) flag.
 //
-//               The example demonstrates user-defined bilinear and nonlinear
-//               form integrators for systems of equations that are defined with
-//               block vectors, and how these are used with an operator for
-//               explicit time integrators. In this case the system also
-//               involves an external approximate Riemann solver for the DG
-//               interface flux. It also demonstrates how to use GLVis for
-//               in-situ visualization of vector grid functions.
+//               The example demonstrates usage of DGHyperbolicConservationLaws
+//               that wraps NonlinearFormIntegrators containing element and face
+//               integration schemes. In this case the system also involves an
+//               external approximate Riemann solver for the DG interface flux.
+//               It also demonstrates how to use GLVis for in-situ visualization
+//               of vector grid function and how to set top-view.
 //
 //               We recommend viewing examples 9, 14 and 17 before viewing this
 //               example.
 
-#include "mfem.hpp"
 #include <fstream>
-#include <sstream>
 #include <iostream>
+#include <sstream>
+#include <cmath>
 
-// Classes FE_Evolution, RiemannSolver, and FaceIntegrator
+#include "mfem.hpp"
+
+// Classes HyperbolicConservationLaws, RiemannSolver, and FaceIntegrator
 // shared between the serial and parallel version of the example.
-#include "ex18.hpp"
+#include "fem/hyperbolic_conservation_laws.hpp"
 
-// Choice for the problem setup. See InitialCondition in ex18.hpp.
-int problem;
+using namespace std;
+using namespace mfem;
 
-// Equation constant parameters.
-const int num_equation = 4;
-const double specific_heat_ratio = 1.4;
-const double gas_constant = 1.0;
+void EulerMesh(const int problem, const char **mesh_file);
 
-// Maximum characteristic speed (updated by integrators)
-double max_char_speed;
+VectorFunctionCoefficient EulerInitialCondition(const int problem,
+                                                const double specific_heat_ratio,
+                                                const double gas_constant);
 
 int main(int argc, char *argv[])
 {
-   // 1. Initialize MPI and HYPRE.
    Mpi::Init(argc, argv);
+   const int numProcs = Mpi::WorldSize();
+   const int myRank = Mpi::WorldRank();
    Hypre::Init();
 
-   // 2. Parse command-line options.
-   problem = 1;
-   const char *mesh_file = "../data/periodic-square.mesh";
+   // 1. Parse command-line options.
+   int problem = 1;
+   const double specific_heat_ratio = 1.4;
+   const double gas_constant = 1.0;
+
+   const char *mesh_file = "";
+   int IntOrderOffset = 3;
    int ser_ref_levels = 0;
-   int par_ref_levels = 1;
+   int par_ref_levels = 2;
    int order = 3;
    int ode_solver_type = 4;
    double t_final = 2.0;
@@ -78,26 +84,22 @@ int main(int argc, char *argv[])
    int vis_steps = 50;
 
    int precision = 8;
-   cout.precision(precision);
+   out.precision(precision);
 
    OptionsParser args(argc, argv);
-   args.AddOption(&mesh_file, "-m", "--mesh",
-                  "Mesh file to use.");
+   args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
    args.AddOption(&problem, "-p", "--problem",
                   "Problem setup to use. See options in velocity_function().");
-   args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
-                  "Number of times to refine the mesh uniformly before parallel"
-                  " partitioning, -1 for auto.");
-   args.AddOption(&par_ref_levels, "-rp", "--refine-parallel",
-                  "Number of times to refine the mesh uniformly after parallel"
-                  " partitioning.");
+   args.AddOption(&ser_ref_levels, "-rs", "--serial-refine",
+                  "Number of times to refine the serial mesh uniformly.");
+   args.AddOption(&par_ref_levels, "-rp", "--parallel-refine",
+                  "Number of times to refine the parallel mesh uniformly.");
    args.AddOption(&order, "-o", "--order",
                   "Order (degree) of the finite elements.");
    args.AddOption(&ode_solver_type, "-s", "--ode-solver",
                   "ODE solver: 1 - Forward Euler,\n\t"
                   "            2 - RK2 SSP, 3 - RK3 SSP, 4 - RK4, 6 - RK6.");
-   args.AddOption(&t_final, "-tf", "--t-final",
-                  "Final time; start time is 0.");
+   args.AddOption(&t_final, "-tf", "--t-final", "Final time; start time is 0.");
    args.AddOption(&dt, "-dt", "--time-step",
                   "Time step. Positive number skips CFL timestep calculation.");
    args.AddOption(&cfl, "-c", "--cfl-number",
@@ -111,55 +113,85 @@ int main(int argc, char *argv[])
    args.Parse();
    if (!args.Good())
    {
-      if (Mpi::Root()) { args.PrintUsage(cout); }
+      if (Mpi::Root()) { args.PrintUsage(out); }
       return 1;
    }
-   if (Mpi::Root()) { args.PrintOptions(cout); }
-
-   // 3. Read the mesh from the given mesh file. This example requires a 2D
-   //    periodic mesh, such as ../data/periodic-square.mesh.
-   Mesh mesh(mesh_file, 1, 1);
-   const int dim = mesh.Dimension();
-
-   MFEM_ASSERT(dim == 2, "Need a two-dimensional mesh for the problem definition");
-
-   // 4. Define the ODE solver used for time integration. Several explicit
-   //    Runge-Kutta methods are available.
-   ODESolver *ode_solver = NULL;
-   switch (ode_solver_type)
+   // When the user does not provide mesh file,
+   // use the default mesh file for the problem.
+   if ((mesh_file == NULL) || (mesh_file[0] == '\0'))    // if NULL or empty
    {
-      case 1: ode_solver = new ForwardEulerSolver; break;
-      case 2: ode_solver = new RK2Solver(1.0); break;
-      case 3: ode_solver = new RK3SSPSolver; break;
-      case 4: ode_solver = new RK4Solver; break;
-      case 6: ode_solver = new RK6Solver; break;
-      default:
-         if (Mpi::Root())
-         {
-            cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
-         }
-         return 3;
+      EulerMesh(problem, &mesh_file);  // get default mesh file name
    }
+   if (Mpi::Root()) { args.PrintOptions(out); }
 
-   // 5. 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.
+   // 2. Read the mesh from the given mesh file.
+   Mesh mesh = Mesh(mesh_file);
+   const int dim = mesh.Dimension();
+   const int num_equations = dim + 2;
+
+   if (problem == 5)
+   {
+      mesh.Transform([](const Vector &x, Vector &y)
+      {
+         y = x;
+         y *= 0.5;
+      });
+   }
+   // perform uniform refine
    for (int lev = 0; lev < ser_ref_levels; lev++)
    {
       mesh.UniformRefinement();
    }
 
-   // 6. Define a parallel mesh by a partitioning of the serial mesh. Refine
-   //    this mesh further in parallel to increase the resolution. Once the
-   //    parallel mesh is defined, the serial mesh can be deleted.
-   ParMesh pmesh(MPI_COMM_WORLD, mesh);
+   if (numProcs > mesh.GetNE())
+   {
+      if (Mpi::Root())
+      {
+         mfem_warning(
+            "The number of processor is larger than the number of elements.\n"
+            "Refine serial meshes until the number of elements is large enough");
+      }
+      while (mesh.GetNE() < numProcs)
+      {
+         mesh.UniformRefinement();
+      }
+   }
+   if (dim > 1) { mesh.EnsureNCMesh(); }
+
+   ParMesh pmesh = ParMesh(MPI_COMM_WORLD, mesh);
    mesh.Clear();
    for (int lev = 0; lev < par_ref_levels; lev++)
    {
       pmesh.UniformRefinement();
    }
+   if (dim > 1) { pmesh.EnsureNCMesh(); }
 
-   // 7. Define the discontinuous DG finite element space of the given
+   // 3. Define the ODE solver used for time integration. Several explicit
+   //    Runge-Kutta methods are available.
+   ODESolver *ode_solver = NULL;
+   switch (ode_solver_type)
+   {
+      case 1:
+         ode_solver = new ForwardEulerSolver;
+         break;
+      case 2:
+         ode_solver = new RK2Solver(1.0);
+         break;
+      case 3:
+         ode_solver = new RK3SSPSolver;
+         break;
+      case 4:
+         ode_solver = new RK4Solver;
+         break;
+      case 6:
+         ode_solver = new RK6Solver;
+         break;
+      default:
+         out << "Unknown ODE solver type: " << ode_solver_type << '\n';
+         return 3;
+   }
+
+   // 4. Define the discontinuous DG finite element space of the given
    //    polynomial order on the refined mesh.
    DG_FECollection fec(order, dim);
    // Finite element space for a scalar (thermodynamic quantity)
@@ -167,105 +199,84 @@ int main(int argc, char *argv[])
    // Finite element space for a mesh-dim vector quantity (momentum)
    ParFiniteElementSpace dfes(&pmesh, &fec, dim, Ordering::byNODES);
    // Finite element space for all variables together (total thermodynamic state)
-   ParFiniteElementSpace vfes(&pmesh, &fec, num_equation, Ordering::byNODES);
+   ParFiniteElementSpace vfes(&pmesh, &fec, num_equations, Ordering::byNODES);
 
    // This example depends on this ordering of the space.
    MFEM_ASSERT(fes.GetOrdering() == Ordering::byNODES, "");
 
-   HYPRE_BigInt glob_size = vfes.GlobalTrueVSize();
    if (Mpi::Root())
    {
-      cout << "Number of unknowns: " << glob_size << endl;
+      out << "Number of unknowns: " << vfes.GetVSize() << endl;
    }
 
-   // 8. Define the initial conditions, save the corresponding mesh and grid
+   // 6. Define the initial conditions, save the corresponding mesh and grid
    //    functions to a file. This can be opened with GLVis with the -gc option.
-
-   // The solution u has components {density, x-momentum, y-momentum, energy}.
-   // These are stored contiguously in the BlockVector u_block.
-   Array<int> offsets(num_equation + 1);
-   for (int k = 0; k <= num_equation; k++) { offsets[k] = k * vfes.GetNDofs(); }
-   BlockVector u_block(offsets);
-
-   // Momentum grid function on dfes for visualization.
-   ParGridFunction mom(&dfes, u_block.GetData() + offsets[1]);
-
    // Initialize the state.
-   VectorFunctionCoefficient u0(num_equation, InitialCondition);
-   ParGridFunction sol(&vfes, u_block.GetData());
+   VectorFunctionCoefficient u0 = EulerInitialCondition(problem,
+                                                        specific_heat_ratio, gas_constant);
+   ParGridFunction sol(&vfes);
    sol.ProjectCoefficient(u0);
 
    // Output the initial solution.
    {
       ostringstream mesh_name;
-      mesh_name << "vortex-mesh." << setfill('0')
-                << setw(6) << Mpi::WorldRank();
+      mesh_name << "euler-mesh." << setfill('0') << setw(6) << Mpi::WorldRank();
       ofstream mesh_ofs(mesh_name.str().c_str());
       mesh_ofs.precision(precision);
       mesh_ofs << pmesh;
 
-      for (int k = 0; k < num_equation; k++)
+      for (int k = 0; k < num_equations; k++)
       {
-         ParGridFunction uk(&fes, u_block.GetBlock(k));
+         ParGridFunction uk(&fes, sol.GetData() + k * fes.GetNDofs());
          ostringstream sol_name;
-         sol_name << "vortex-" << k << "-init."
-                  << setfill('0') << setw(6) << Mpi::WorldRank();
+         sol_name << "euler-" << k << "-init." << setfill('0') << setw(6)
+                  << Mpi::WorldRank();
          ofstream sol_ofs(sol_name.str().c_str());
          sol_ofs.precision(precision);
          sol_ofs << uk;
       }
    }
 
-   // 9. Set up the nonlinear form corresponding to the DG discretization of the
+   // 7. Set up the nonlinear form corresponding to the DG discretization of the
    //    flux divergence, and assemble the corresponding mass matrix.
-   MixedBilinearForm Aflux(&dfes, &fes);
-   Aflux.AddDomainIntegrator(new TransposeIntegrator(new GradientIntegrator()));
-   Aflux.Assemble();
-
-   ParNonlinearForm A(&vfes);
-   RiemannSolver rsolver;
-   A.AddInteriorFaceIntegrator(new FaceIntegrator(rsolver, dim));
-
-   // 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 euler(vfes, A, Aflux.SpMat());
+   RiemannSolver *numericalFlux = new RusanovFlux();
+   DGHyperbolicConservationLaws euler = getEulerSystem(
+                                           &vfes, numericalFlux, specific_heat_ratio, IntOrderOffset);
 
    // Visualize the density
    socketstream sout;
    if (visualization)
    {
       char vishost[] = "localhost";
-      int  visport   = 19916;
+      int visport = 19916;
 
-      MPI_Barrier(pmesh.GetComm());
       sout.open(vishost, visport);
       if (!sout)
       {
-         if (Mpi::Root())
-         {
-            cout << "Unable to connect to GLVis server at "
-                 << vishost << ':' << visport << endl;
-         }
          visualization = false;
          if (Mpi::Root())
          {
-            cout << "GLVis visualization disabled.\n";
+            out << "Unable to connect to GLVis server at " << vishost << ':'
+                << visport << endl;
+            out << "GLVis visualization disabled.\n";
          }
       }
       else
       {
-         sout << "parallel " << Mpi::WorldSize()
-              << " " << Mpi::WorldRank() << "\n";
+         ParGridFunction mom(&dfes, sol.GetData());
+         sout << "parallel " << numProcs << " " << myRank << "\n";
          sout.precision(precision);
          sout << "solution\n" << pmesh << mom;
+         sout << "view 0 0\n";  // view from top
+         sout << "keys jlm\n";  // turn off perspective and light
          sout << "pause\n";
          sout << flush;
          if (Mpi::Root())
          {
-            cout << "GLVis visualization paused."
-                 << " Press space (in the GLVis window) to resume it.\n";
+            out << "GLVis visualization paused."
+                << " Press space (in the GLVis window) to resume it.\n";
          }
+         MPI_Barrier(pmesh.GetComm());
       }
    }
 
@@ -273,12 +284,11 @@ int main(int argc, char *argv[])
    double hmin;
    if (cfl > 0)
    {
-      double my_hmin = pmesh.GetElementSize(0, 1);
+      double my_hmin = pmesh.GetNE() > 0 ? pmesh.GetElementSize(0, 1) : INFINITY;
       for (int i = 1; i < pmesh.GetNE(); i++)
       {
          my_hmin = min(pmesh.GetElementSize(i, 1), my_hmin);
       }
-      // Reduce to find the global minimum element size
       MPI_Allreduce(&my_hmin, &hmin, 1, MPI_DOUBLE, MPI_MIN, pmesh.GetComm());
    }
 
@@ -294,52 +304,47 @@ int main(int argc, char *argv[])
    {
       // Find a safe dt, using a temporary vector. Calling Mult() computes the
       // maximum char speed at all quadrature points on all faces.
-      max_char_speed = 0.;
       Vector z(sol.Size());
-      A.Mult(sol, z);
-      // Reduce to find the global maximum wave speed
-      {
-         double all_max_char_speed;
-         MPI_Allreduce(&max_char_speed, &all_max_char_speed,
-                       1, MPI_DOUBLE, MPI_MAX, pmesh.GetComm());
-         max_char_speed = all_max_char_speed;
-      }
-      dt = cfl * hmin / max_char_speed / (2*order+1);
+      euler.Mult(sol, z);
+
+      double max_char_speed;
+      double my_max_char_speed = euler.getMaxCharSpeed();
+      MPI_Allreduce(&my_max_char_speed, &max_char_speed, 1, MPI_DOUBLE, MPI_MAX,
+                    pmesh.GetComm());
+      dt = cfl * hmin / max_char_speed / (2 * order + 1);
    }
 
    // Integrate in time.
    bool done = false;
-   for (int ti = 0; !done; )
+   for (int ti = 0; !done;)
    {
       double dt_real = min(dt, t_final - t);
 
       ode_solver->Step(sol, t, dt_real);
       if (cfl > 0)
       {
-         // Reduce to find the global maximum wave speed
-         {
-            double all_max_char_speed;
-            MPI_Allreduce(&max_char_speed, &all_max_char_speed,
-                          1, MPI_DOUBLE, MPI_MAX, pmesh.GetComm());
-            max_char_speed = all_max_char_speed;
-         }
-         dt = cfl * hmin / max_char_speed / (2*order+1);
+         double max_char_speed;
+         double my_max_char_speed = euler.getMaxCharSpeed();
+         MPI_Allreduce(&my_max_char_speed, &max_char_speed, 1, MPI_DOUBLE, MPI_MAX,
+                       pmesh.GetComm());
+         dt = cfl * hmin / max_char_speed / (2 * order + 1);
       }
       ti++;
 
-      done = (t >= t_final - 1e-8*dt);
+      done = (t >= t_final - 1e-8 * dt);
       if (done || ti % vis_steps == 0)
       {
          if (Mpi::Root())
          {
-            cout << "time step: " << ti << ", time: " << t << endl;
+            out << "time step: " << ti << ", time: " << t << endl;
          }
          if (visualization)
          {
-            MPI_Barrier(pmesh.GetComm());
-            sout << "parallel " << Mpi::WorldSize()
-                 << " " << Mpi::WorldRank() << "\n";
+            ParGridFunction mom(&dfes, sol.GetData());
+            sout << "parallel " << numProcs << " " << myRank << "\n";
             sout << "solution\n" << pmesh << mom << flush;
+            sout << "window_title 't = " << t << "'";
+            MPI_Barrier(pmesh.GetComm());
          }
       }
    }
@@ -347,34 +352,229 @@ int main(int argc, char *argv[])
    tic_toc.Stop();
    if (Mpi::Root())
    {
-      cout << " done, " << tic_toc.RealTime() << "s." << endl;
+      out << " done, " << tic_toc.RealTime() << "s." << endl;
    }
 
-   // 11. Save the final solution. This output can be viewed later using GLVis:
-   //     "glvis -np 4 -m vortex-mesh -g vortex-1-final".
-   for (int k = 0; k < num_equation; k++)
+   // 9. Save the final solution. This output can be viewed later using GLVis:
+   //    "glvis -m euler.mesh -g euler-1-final.gf".
    {
-      ParGridFunction uk(&fes, u_block.GetBlock(k));
-      ostringstream sol_name;
-      sol_name << "vortex-" << k << "-final."
-               << setfill('0') << setw(6) << Mpi::WorldRank();
-      ofstream sol_ofs(sol_name.str().c_str());
-      sol_ofs.precision(precision);
-      sol_ofs << uk;
-   }
+      ostringstream mesh_name;
+      mesh_name << "euler-mesh-final." << setfill('0') << setw(6)
+                << Mpi::WorldRank();
+      ofstream mesh_ofs(mesh_name.str().c_str());
+      mesh_ofs.precision(precision);
+      mesh_ofs << pmesh;
 
-   // 12. Compute the L2 solution error summed for all components.
-   if (t_final == 2.0)
-   {
-      const double error = sol.ComputeLpError(2, u0);
-      if (Mpi::Root())
+      for (int k = 0; k < num_equations; k++)
       {
-         cout << "Solution error: " << error << endl;
+         ParGridFunction uk(&fes, sol.GetData() + k * fes.GetNDofs());
+         ostringstream sol_name;
+         sol_name << "euler-" << k << "-final." << setfill('0') << setw(6)
+                  << Mpi::WorldRank();
+         ofstream sol_ofs(sol_name.str().c_str());
+         sol_ofs.precision(precision);
+         sol_ofs << uk;
       }
    }
 
+   // 10. Compute the L2 solution error summed for all components.
+   //   if (t_final == 2.0) {
+   const double error = sol.ComputeLpError(2, u0);
+   if (Mpi::Root())
+   {
+      out << "Solution error: " << error << endl;
+   }
+
    // Free the used memory.
    delete ode_solver;
 
    return 0;
 }
+
+void EulerMesh(const int problem, const char **mesh_file)
+{
+   switch (problem)
+   {
+      case 1:
+         *mesh_file = "../data/periodic-square-4x4.mesh";
+         break;
+      case 2:
+         *mesh_file = "../data/periodic-square-4x4.mesh";
+         break;
+      case 3:
+         *mesh_file = "../data/periodic-square-4x4.mesh";
+         break;
+      case 4:
+         *mesh_file = "../data/periodic-segment.mesh";
+         break;
+      case 5:
+         *mesh_file = "../data/periodic-square-4x4.mesh";
+         break;
+      default:
+         throw invalid_argument("Default mesh is undefined");
+   }
+}
+
+// Initial condition
+VectorFunctionCoefficient EulerInitialCondition(const int problem,
+                                                const double specific_heat_ratio,
+                                                const double gas_constant)
+{
+   switch (problem)
+   {
+      case 1: // fast moving vortex
+         return VectorFunctionCoefficient(4, [specific_heat_ratio,
+                                              gas_constant](const Vector &x, Vector &y)
+         {
+            MFEM_ASSERT(x.Size() == 2, "");
+
+            double radius = 0, Minf = 0, beta = 0;
+            // "Fast euler"
+            radius = 0.2;
+            Minf = 0.5;
+            beta = 1. / 5.;
+
+            const double xc = 0.0, yc = 0.0;
+
+            // Nice units
+            const double vel_inf = 1.;
+            const double den_inf = 1.;
+
+            // Derive remainder of background state from this and Minf
+            const double pres_inf = (den_inf / specific_heat_ratio) *
+                                    (vel_inf / Minf) * (vel_inf / Minf);
+            const double temp_inf = pres_inf / (den_inf * gas_constant);
+
+            double r2rad = 0.0;
+            r2rad += (x(0) - xc) * (x(0) - xc);
+            r2rad += (x(1) - yc) * (x(1) - yc);
+            r2rad /= (radius * radius);
+
+            const double shrinv1 = 1.0 / (specific_heat_ratio - 1.);
+
+            const double velX =
+               vel_inf * (1 - beta * (x(1) - yc) / radius * exp(-0.5 * r2rad));
+            const double velY =
+               vel_inf * beta * (x(0) - xc) / radius * exp(-0.5 * r2rad);
+            const double vel2 = velX * velX + velY * velY;
+
+            const double specific_heat =
+               gas_constant * specific_heat_ratio * shrinv1;
+            const double temp = temp_inf - 0.5 * (vel_inf * beta) *
+                                (vel_inf * beta) / specific_heat *
+                                exp(-r2rad);
+
+            const double den = den_inf * pow(temp / temp_inf, shrinv1);
+            const double pres = den * gas_constant * temp;
+            const double energy = shrinv1 * pres / den + 0.5 * vel2;
+
+            y(0) = den;
+            y(1) = den * velX;
+            y(2) = den * velY;
+            y(3) = den * energy;
+         });
+      case 2: // slow moving vortex
+         return VectorFunctionCoefficient(4, [specific_heat_ratio,
+                                              gas_constant](const Vector &x, Vector &y)
+         {
+            MFEM_ASSERT(x.Size() == 2, "");
+
+            double radius = 0, Minf = 0, beta = 0;
+            // "Slow euler"
+            radius = 0.2;
+            Minf = 0.05;
+            beta = 1. / 50.;
+
+            const double xc = 0.0, yc = 0.0;
+
+            // Nice units
+            const double vel_inf = 1.;
+            const double den_inf = 1.;
+
+            // Derive remainder of background state from this and Minf
+            const double pres_inf = (den_inf / specific_heat_ratio) *
+                                    (vel_inf / Minf) * (vel_inf / Minf);
+            const double temp_inf = pres_inf / (den_inf * gas_constant);
+
+            double r2rad = 0.0;
+            r2rad += (x(0) - xc) * (x(0) - xc);
+            r2rad += (x(1) - yc) * (x(1) - yc);
+            r2rad /= (radius * radius);
+
+            const double shrinv1 = 1.0 / (specific_heat_ratio - 1.);
+
+            const double velX =
+               vel_inf * (1 - beta * (x(1) - yc) / radius * exp(-0.5 * r2rad));
+            const double velY =
+               vel_inf * beta * (x(0) - xc) / radius * exp(-0.5 * r2rad);
+            const double vel2 = velX * velX + velY * velY;
+
+            const double specific_heat =
+               gas_constant * specific_heat_ratio * shrinv1;
+            const double temp = temp_inf - 0.5 * (vel_inf * beta) *
+                                (vel_inf * beta) / specific_heat *
+                                exp(-r2rad);
+
+            const double den = den_inf * pow(temp / temp_inf, shrinv1);
+            const double pres = den * gas_constant * temp;
+            const double energy = shrinv1 * pres / den + 0.5 * vel2;
+
+            y(0) = den;
+            y(1) = den * velX;
+            y(2) = den * velY;
+            y(3) = den * energy;
+         });
+      case 3: // moving sine wave
+         return VectorFunctionCoefficient(4, [](const Vector &x, Vector &y)
+         {
+            MFEM_ASSERT(x.Size() == 2, "");
+            const double density = 1.0 + 0.2 * sin(M_PI*(x(0) + x(1)));
+            const double velocity_x = 0.7;
+            const double velocity_y = 0.3;
+            const double pressure = 1.0;
+            const double energy =
+               pressure / (1.4 - 1.0) +
+               density * 0.5 * (velocity_x * velocity_x + velocity_y * velocity_y);
+
+            y(0) = density;
+            y(1) = density * velocity_x;
+            y(2) = density * velocity_y;
+            y(3) = energy;
+         });
+      case 4:
+         return VectorFunctionCoefficient(3, [](const Vector &x, Vector &y)
+         {
+            MFEM_ASSERT(x.Size() == 1, "");
+            const double density = 1.0 + 0.2 * sin(M_PI * 2 * x(0));
+            const double velocity_x = 1.0;
+            const double pressure = 1.0;
+            const double energy =
+               pressure / (1.4 - 1.0) + density * 0.5 * (velocity_x * velocity_x);
+
+            y(0) = density;
+            y(1) = density * velocity_x;
+            y(2) = energy;
+         });
+      case 5:
+         return VectorFunctionCoefficient(4, [](const Vector &x, Vector &y)
+         {
+            MFEM_ASSERT(x.Size() == 2, "");
+            const double L = 1.0;
+            const double density = abs(x(1)) < 0.25 ? 2 : 1;
+            const double velocity_x = abs(x(1)) < 0.25 ? -0.5 : 0.5;
+            const double velocity_y = abs(x(1)) < 0.25 ? 0.01 * sin(M_PI*x(0) / L)
+                                      : 0.01 * sin(M_PI*x(0) / L);
+            const double pressure = abs(x(1)) < 0.25 ? 2.5 : 2.5;
+            const double energy =
+               pressure / (1.4 - 1.0) +
+               density * 0.5 * (velocity_x * velocity_x + velocity_y * velocity_y);
+
+            y(0) = density;
+            y(1) = density * velocity_x;
+            y(2) = density * velocity_y;
+            y(3) = energy;
+         });
+      default:
+         throw invalid_argument("Problem Undefined");
+   }
+}

examples/makefile

@@ -90,11 +90,10 @@ $(SUBDIRS_TPRINT):
 	@$(MAKE) -C $(@D) $(@F)
 
 # Additional dependencies
-ex18: $(SRC)ex18.hpp
 ex33: $(SRC)ex33.hpp
+ex34: $(SRC)ex34.hpp
 
 ifeq ($(MFEM_USE_MPI),YES)
-ex18p: $(SRC)ex18.hpp
 ex33p: $(SRC)ex33.hpp
 endif
 

examples/proxGalerkinAdv.cpp

@@ -0,0 +1,723 @@
+//                                MFEM Example 9
+//
+// Compile with: make ex9
+//
+// Sample runs:
+//    ex9 -m ../data/periodic-segment.mesh -p 0 -r 2 -dt 0.005
+//    ex9 -m ../data/periodic-square.mesh -p 0 -r 2 -dt 0.01 -tf 10
+//    ex9 -m ../data/periodic-hexagon.mesh -p 0 -r 2 -dt 0.01 -tf 10
+//    ex9 -m ../data/periodic-square.mesh -p 1 -r 2 -dt 0.005 -tf 9
+//    ex9 -m ../data/periodic-hexagon.mesh -p 1 -r 2 -dt 0.005 -tf 9
+//    ex9 -m ../data/amr-quad.mesh -p 1 -r 2 -dt 0.002 -tf 9
+//    ex9 -m ../data/amr-quad.mesh -p 1 -r 2 -dt 0.02 -s 13 -tf 9
+//    ex9 -m ../data/star-q3.mesh -p 1 -r 2 -dt 0.005 -tf 9
+//    ex9 -m ../data/star-mixed.mesh -p 1 -r 2 -dt 0.005 -tf 9
+//    ex9 -m ../data/disc-nurbs.mesh -p 1 -r 3 -dt 0.005 -tf 9
+//    ex9 -m ../data/disc-nurbs.mesh -p 2 -r 3 -dt 0.005 -tf 9
+//    ex9 -m ../data/periodic-square.mesh -p 3 -r 4 -dt 0.0025 -tf 9 -vs 20
+//    ex9 -m ../data/periodic-cube.mesh -p 0 -r 2 -o 2 -dt 0.02 -tf 8
+//    ex9 -m ../data/periodic-square.msh -p 0 -r 2 -dt 0.005 -tf 2
+//    ex9 -m ../data/periodic-cube.msh -p 0 -r 1 -o 2 -tf 2
+//
+// Device sample runs:
+//    ex9 -pa
+//    ex9 -ea
+//    ex9 -fa
+//    ex9 -pa -m ../data/periodic-cube.mesh
+//    ex9 -pa -m ../data/periodic-cube.mesh -d cuda
+//    ex9 -ea -m ../data/periodic-cube.mesh -d cuda
+//    ex9 -fa -m ../data/periodic-cube.mesh -d cuda
+//    ex9 -pa -m ../data/amr-quad.mesh -p 1 -r 2 -dt 0.002 -tf 9 -d cuda
+//
+// Description:  This example code solves the time-dependent advection equation
+//               du/dt + v.grad(u) = 0, where v is a given fluid velocity, and
+//               u0(x)=u(0,x) is a given initial condition.
+//
+//               The example demonstrates the use of Discontinuous Galerkin (DG)
+//               bilinear forms in MFEM (face integrators), the use of implicit
+//               and explicit ODE time integrators, the definition of periodic
+//               boundary conditions through periodic meshes, as well as the use
+//               of GLVis for persistent visualization of a time-evolving
+//               solution. The saving of time-dependent data files for external
+//               visualization with VisIt (visit.llnl.gov) and ParaView
+//               (paraview.org) is also illustrated.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+#include <algorithm>
+
+using namespace std;
+using namespace mfem;
+
+// Choice for the problem setup. The fluid velocity, initial condition and
+// inflow boundary condition are chosen based on this parameter.
+int problem;
+
+// Velocity coefficient
+void velocity_function(const Vector &x, Vector &v);
+
+// Initial condition
+double u0_function(const Vector &x);
+
+// Inflow boundary condition
+double inflow_function(const Vector &x);
+
+// Mesh bounding box
+Vector bb_min, bb_max;
+
+
+class DG_Solver : public Solver
+{
+private:
+   SparseMatrix &M, &K, A;
+   GMRESSolver linear_solver;
+   BlockILU prec;
+   double dt;
+public:
+   DG_Solver(SparseMatrix &M_, SparseMatrix &K_, const FiniteElementSpace &fes)
+      : M(M_),
+        K(K_),
+        prec(fes.GetFE(0)->GetDof(),
+             BlockILU::Reordering::MINIMUM_DISCARDED_FILL),
+        dt(-1.0)
+   {
+      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);
+   }
+
+   void SetTimeStep(double dt_)
+   {
+      if (dt_ != dt)
+      {
+         dt = dt_;
+         // Form operator A = M - dt*K
+         A = K;
+         A *= -dt;
+         A += M;
+
+         // this will also call SetOperator on the preconditioner
+         linear_solver.SetOperator(A);
+      }
+   }
+
+   void SetOperator(const Operator &op)
+   {
+      linear_solver.SetOperator(op);
+   }
+
+   virtual void Mult(const Vector &x, Vector &y) const
+   {
+      linear_solver.Mult(x, y);
+   }
+};
+
+/** 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
+    be written as a general ODE, du/dt = M^{-1} (K u + b), and this class is
+    used to evaluate the right-hand side. */
+class FE_Evolution : public TimeDependentOperator
+{
+private:
+   BilinearForm &M, &K;
+   const Vector &b;
+   Solver *M_prec;
+   CGSolver M_solver;
+   DG_Solver *dg_solver;
+   FiniteElementSpace *fes;
+
+   GridFunction latent;
+   GridFunction latent_k;
+   GridFunction xnew;
+
+   GridFunctionCoefficient latent_cf;
+   GridFunctionCoefficient latent_k_cf;
+   GridFunctionCoefficient xnew_cf;
+   Vector delta_k, delta_latent;
+
+   TransformedCoefficient expLatent_cf;
+   SumCoefficient diff_latent_cf;
+   TransformedCoefficient expResidual_cf;
+
+   LinearForm delta_latent_form;
+   NonlinearForm invExpLatentForm;
+   LinearForm newtonRHS;
+
+   mutable Vector z;
+
+public:
+   FE_Evolution(BilinearForm &M_, BilinearForm &K_, const Vector &b_);
+
+   virtual void Mult(const Vector &x, Vector &y) const;
+   virtual void ImplicitSolve(const double dt, const Vector &x, Vector &k);
+
+   virtual ~FE_Evolution();
+   bool postprocess = false;
+};
+
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   problem = 1;
+   const char *mesh_file = "../data/periodic-square-4x4.mesh";
+   int ref_levels = 4;
+   int order = 1;
+   bool pa = false;
+   bool ea = false;
+   bool fa = false;
+   const char *device_config = "cpu";
+   int ode_solver_type = 11;
+   double t_final = 10.0;
+   double dt = 0.01;
+   bool visualization = true;
+   bool visit = false;
+   bool paraview = false;
+   bool binary = false;
+   int vis_steps = 5;
+   bool applyPostprocess = false;
+
+   int precision = 8;
+   cout.precision(precision);
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&problem, "-p", "--problem",
+                  "Problem setup to use. See options in velocity_function().");
+   args.AddOption(&ref_levels, "-r", "--refine",
+                  "Number of times to refine the mesh uniformly.");
+   args.AddOption(&order, "-o", "--order",
+                  "Order (degree) of the finite elements.");
+   args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
+                  "--no-partial-assembly", "Enable Partial Assembly.");
+   args.AddOption(&ea, "-ea", "--element-assembly", "-no-ea",
+                  "--no-element-assembly", "Enable Element Assembly.");
+   args.AddOption(&fa, "-fa", "--full-assembly", "-no-fa",
+                  "--no-full-assembly", "Enable Full Assembly.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&ode_solver_type, "-s", "--ode-solver",
+                  "ODE solver: 1 - Forward Euler,\n\t"
+                  "            2 - RK2 SSP, 3 - RK3 SSP, 4 - RK4, 6 - RK6,\n\t"
+                  "            11 - Backward Euler,\n\t"
+                  "            12 - SDIRK23 (L-stable), 13 - SDIRK33,\n\t"
+                  "            22 - Implicit Midpoint Method,\n\t"
+                  "            23 - SDIRK23 (A-stable), 24 - SDIRK34");
+   args.AddOption(&t_final, "-tf", "--t-final",
+                  "Final time; start time is 0.");
+   args.AddOption(&dt, "-dt", "--time-step",
+                  "Time step.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&visit, "-visit", "--visit-datafiles", "-no-visit",
+                  "--no-visit-datafiles",
+                  "Save data files for VisIt (visit.llnl.gov) visualization.");
+   args.AddOption(&paraview, "-paraview", "--paraview-datafiles", "-no-paraview",
+                  "--no-paraview-datafiles",
+                  "Save data files for ParaView (paraview.org) visualization.");
+   args.AddOption(&binary, "-binary", "--binary-datafiles", "-ascii",
+                  "--ascii-datafiles",
+                  "Use binary (Sidre) or ascii format for VisIt data files.");
+   args.AddOption(&vis_steps, "-vs", "--visualization-steps",
+                  "Visualize every n-th timestep.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   Device device(device_config);
+   device.Print();
+
+   // 2. Read the mesh from the given mesh file. We can handle geometrically
+   //    periodic meshes in this code.
+   Mesh mesh(mesh_file, 1, 1);
+   int dim = mesh.Dimension();
+
+   // 3. Define the ODE solver used for time integration. Several explicit
+   //    Runge-Kutta methods are available.
+   ODESolver *ode_solver = NULL;
+   switch (ode_solver_type)
+   {
+      // Explicit methods
+      case 1: ode_solver = new ForwardEulerSolver; break;
+      case 2: ode_solver = new RK2Solver(1.0); break;
+      case 3: ode_solver = new RK3SSPSolver; break;
+      case 4: ode_solver = new RK4Solver; break;
+      case 6: ode_solver = new RK6Solver; break;
+      // Implicit (L-stable) methods
+      case 11: ode_solver = new BackwardEulerSolver; break;
+      case 12: ode_solver = new SDIRK23Solver(2); break;
+      case 13: ode_solver = new SDIRK33Solver; break;
+      // Implicit A-stable methods (not L-stable)
+      case 22: ode_solver = new ImplicitMidpointSolver; break;
+      case 23: ode_solver = new SDIRK23Solver; break;
+      case 24: ode_solver = new SDIRK34Solver; break;
+
+      default:
+         cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
+         return 3;
+   }
+
+   // 4. Refine the mesh to increase the resolution. In this example we do
+   //    'ref_levels' of uniform refinement, where 'ref_levels' is a
+   //    command-line parameter. If the mesh is of NURBS type, we convert it to
+   //    a (piecewise-polynomial) high-order mesh.
+   for (int lev = 0; lev < ref_levels; lev++)
+   {
+      mesh.UniformRefinement();
+   }
+   if (mesh.NURBSext)
+   {
+      mesh.SetCurvature(max(order, 1));
+   }
+   mesh.GetBoundingBox(bb_min, bb_max, max(order, 1));
+
+   // 5. Define the discontinuous DG finite element space of the given
+   //    polynomial order on the refined mesh.
+   DG_FECollection fec(order, dim, BasisType::GaussLobatto);
+   FiniteElementSpace fes(&mesh, &fec);
+
+   cout << "Number of unknowns: " << fes.GetVSize() << endl;
+
+   // 6. Set up and assemble the bilinear and linear forms corresponding to the
+   //    DG discretization. The DGTraceIntegrator involves integrals over mesh
+   //    interior faces.
+   VectorFunctionCoefficient velocity(dim, velocity_function);
+   FunctionCoefficient inflow(inflow_function);
+   FunctionCoefficient u0(u0_function);
+
+   BilinearForm m(&fes);
+   BilinearForm k(&fes);
+   if (pa)
+   {
+      m.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+      k.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+   }
+   else if (ea)
+   {
+      m.SetAssemblyLevel(AssemblyLevel::ELEMENT);
+      k.SetAssemblyLevel(AssemblyLevel::ELEMENT);
+   }
+   else if (fa)
+   {
+      m.SetAssemblyLevel(AssemblyLevel::FULL);
+      k.SetAssemblyLevel(AssemblyLevel::FULL);
+   }
+   m.AddDomainIntegrator(new MassIntegrator);
+   constexpr double alpha = -1.0;
+   k.AddDomainIntegrator(new ConvectionIntegrator(velocity, alpha));
+   k.AddInteriorFaceIntegrator(
+      new NonconservativeDGTraceIntegrator(velocity, alpha));
+   k.AddBdrFaceIntegrator(
+      new NonconservativeDGTraceIntegrator(velocity, alpha));
+
+   LinearForm b(&fes);
+   b.AddBdrFaceIntegrator(
+      new BoundaryFlowIntegrator(inflow, velocity, alpha));
+
+   m.Assemble();
+   int skip_zeros = 0;
+   k.Assemble(skip_zeros);
+   b.Assemble();
+   m.Finalize();
+   k.Finalize(skip_zeros);
+
+   // 7. Define the initial conditions, save the corresponding grid function to
+   //    a file and (optionally) save data in the VisIt format and initialize
+   //    GLVis visualization.
+   GridFunction u(&fes);
+   u.ProjectCoefficient(u0);
+
+   {
+      ofstream omesh("ex9.mesh");
+      omesh.precision(precision);
+      mesh.Print(omesh);
+      ofstream osol("ex9-init.gf");
+      osol.precision(precision);
+      u.Save(osol);
+   }
+
+   // Create data collection for solution output: either VisItDataCollection for
+   // ascii data files, or SidreDataCollection for binary data files.
+   DataCollection *dc = NULL;
+   if (visit)
+   {
+      if (binary)
+      {
+#ifdef MFEM_USE_SIDRE
+         dc = new SidreDataCollection("Example9", &mesh);
+#else
+         MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
+#endif
+      }
+      else
+      {
+         dc = new VisItDataCollection("Example9", &mesh);
+         dc->SetPrecision(precision);
+      }
+      dc->RegisterField("solution", &u);
+      dc->SetCycle(0);
+      dc->SetTime(0.0);
+      dc->Save();
+   }
+
+   ParaViewDataCollection *pd = NULL;
+   if (paraview)
+   {
+      pd = new ParaViewDataCollection("Example9", &mesh);
+      pd->SetPrefixPath("ParaView");
+      pd->RegisterField("solution", &u);
+      pd->SetLevelsOfDetail(order);
+      pd->SetDataFormat(VTKFormat::BINARY);
+      pd->SetHighOrderOutput(true);
+      pd->SetCycle(0);
+      pd->SetTime(0.0);
+      pd->Save();
+   }
+
+   socketstream sout;
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      sout.open(vishost, visport);
+      if (!sout)
+      {
+         cout << "Unable to connect to GLVis server at "
+              << vishost << ':' << visport << endl;
+         visualization = false;
+         cout << "GLVis visualization disabled.\n";
+      }
+      else
+      {
+         sout.precision(precision);
+         sout << "solution\n" << mesh << u;
+         sout << "pause\n";
+         sout << flush;
+         cout << "GLVis visualization paused."
+              << " Press space (in the GLVis window) to resume it.\n";
+      }
+   }
+
+   // 8. 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);
+   adv.postprocess = applyPostprocess;
+
+   double t = 0.0;
+   adv.SetTime(t);
+   ode_solver->Init(adv);
+
+   bool done = false;
+   for (int ti = 0; !done; )
+   {
+      double dt_real = min(dt, t_final - t);
+      ode_solver->Step(u, t, dt_real);
+      ti++;
+
+      done = (t >= t_final - 1e-8*dt);
+
+      if (done || ti % vis_steps == 0)
+      {
+         cout << "time step: " << ti << ", time: " << t << endl;
+
+         if (visualization)
+         {
+            sout << "solution\n" << mesh << u << flush;
+         }
+
+         if (visit)
+         {
+            dc->SetCycle(ti);
+            dc->SetTime(t);
+            dc->Save();
+         }
+
+         if (paraview)
+         {
+            pd->SetCycle(ti);
+            pd->SetTime(t);
+            pd->Save();
+         }
+      }
+   }
+
+   // 9. Save the final solution. This output can be viewed later using GLVis:
+   //    "glvis -m ex9.mesh -g ex9-final.gf".
+   {
+      ofstream osol("ex9-final.gf");
+      osol.precision(precision);
+      u.Save(osol);
+   }
+
+   // 10. Free the used memory.
+   delete ode_solver;
+   delete pd;
+   delete dc;
+
+   return 0;
+}
+
+
+// Implementation of class FE_Evolution
+FE_Evolution::FE_Evolution(BilinearForm &M_, BilinearForm &K_, const Vector &b_)
+   : TimeDependentOperator(M_.Height()), M(M_), K(K_), b(b_), fes(M.FESpace()),
+     latent(fes), latent_k(fes), xnew(fes), latent_cf(&latent),
+     latent_k_cf(&latent_k), xnew_cf(&xnew), delta_k(xnew.Size()),
+     delta_latent(latent.Size()), expLatent_cf(&latent_cf, exp),
+     diff_latent_cf(latent_cf, latent_k_cf, -1.0, 1.0), expResidual_cf(&xnew_cf,
+                                                                       &latent_cf, [](double x, double psi) {return x - exp(psi)*(1.0 - psi);}),
+delta_latent_form(fes), invExpLatentForm(fes), newtonRHS(fes), z(M_.Height())
+{
+
+   delta_latent_form.AddDomainIntegrator(new DomainLFIntegrator(diff_latent_cf));
+   // since DG is local space, we can perform inverse element-wise.
+   invExpLatentForm.AddDomainIntegrator(new InverseIntegrator(new MassIntegrator(
+                                                                 expLatent_cf)));
+   newtonRHS.AddDomainIntegrator(new DomainLFIntegrator(expResidual_cf));
+   Array<int> ess_tdof_list;
+   if (M.GetAssemblyLevel() == AssemblyLevel::LEGACY)
+   {
+      M_prec = new DSmoother(M.SpMat());
+      M_solver.SetOperator(M.SpMat());
+      dg_solver = new DG_Solver(M.SpMat(), K.SpMat(), *M.FESpace());
+   }
+   else
+   {
+      M_prec = new OperatorJacobiSmoother(M, ess_tdof_list);
+      M_solver.SetOperator(M);
+      dg_solver = NULL;
+   }
+   M_solver.SetPreconditioner(*M_prec);
+   M_solver.iterative_mode = false;
+   M_solver.SetRelTol(1e-9);
+   M_solver.SetAbsTol(0.0);
+   M_solver.SetMaxIter(100);
+   M_solver.SetPrintLevel(0);
+}
+
+void FE_Evolution::Mult(const Vector &x, Vector &y) const
+{
+   // y = M^{-1} (K x + b)
+   K.Mult(x, z);
+   z += b;
+   M_solver.Mult(z, y);
+}
+
+void FE_Evolution::ImplicitSolve(const double dt, const Vector &x, Vector &k)
+{
+   MFEM_VERIFY(dg_solver != NULL,
+               "Implicit time integration is not supported with partial assembly");
+
+   k.SetSize(x.Size());
+   k = 0.0;
+
+   latent = x;
+   latent.ApplyMap([](double x) {return x > 0 ? log(x) : -20.0; });
+
+   GridFunction zero_gf(fes);
+   zero_gf = 0.0;
+
+   for (int i=1; i<1e04; i++)
+   {
+      latent_k = latent;
+      for (int j=0; j<1e04; j++)
+      {
+         delta_k = k;
+
+         K.Mult(x, z);
+         z += b;
+         delta_latent_form.Assemble();
+         z.Add(dt*i, delta_latent_form);
+         dg_solver->SetTimeStep(dt);
+         dg_solver->Mult(z, k);
+
+         delta_k -= k;
+
+         xnew = x;
+         xnew.Add(dt, k);
+         latent = x;
+         latent.Add(dt, k);
+         latent.ApplyMap([](double x) {return x > 0 ? log(x) : -20.0;});
+         bool converged = false;
+         for (int jj=0; jj<1e04; jj++)
+         {
+            delta_latent = latent;
+
+            newtonRHS.Assemble();
+            invExpLatentForm.Mult(newtonRHS, latent);
+
+            delta_latent -= latent;
+            if (delta_latent.Normlinf() < 1e-07)
+            {
+               converged = true;
+               break;
+            }
+         }
+         if (!converged) {mfem_warning("Latent failed to converge."); }
+
+         // out << delta_latent.Normlinf() << ", " << delta_k.Normlinf() << std::endl;
+
+         if (delta_latent.Normlinf() < 1e-07 & delta_k.Normlinf() < 1e-07)
+         {
+            out << "\t\tInnerloop convereged in " << j << "steps."<< std::endl;
+            break;
+         }
+      }
+      latent_k -= latent;
+      TransformedCoefficient latent_diff(&latent_k_cf, &latent_cf, [](double x,
+      double xnew) {return exp(x)*(x-xnew);});
+      if (zero_gf.ComputeLpError(infinity(), latent_diff) < 1e-07)
+      {
+         out << "\tOuterloop convereged in " << i << "steps."<< std::endl;
+         break;
+      }
+   }
+   if (postprocess)
+   {
+      k = latent;
+      k.ApplyMap([](double x) {return exp(x); });
+      k -= x;
+      k *= dt;
+   }
+}
+
+FE_Evolution::~FE_Evolution()
+{
+   delete M_prec;
+   delete dg_solver;
+}
+
+// Velocity coefficient
+void velocity_function(const Vector &x, Vector &v)
+{
+   int dim = x.Size();
+
+   // map to the reference [-1,1] domain
+   Vector X(dim);
+   for (int i = 0; i < dim; i++)
+   {
+      double center = (bb_min[i] + bb_max[i]) * 0.5;
+      X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
+   }
+
+   switch (problem)
+   {
+      case 0:
+      {
+         // Translations in 1D, 2D, and 3D
+         switch (dim)
+         {
+            case 1: v(0) = 1.0; break;
+            case 2: v(0) = sqrt(2./3.); v(1) = sqrt(1./3.); break;
+            case 3: v(0) = sqrt(3./6.); v(1) = sqrt(2./6.); v(2) = sqrt(1./6.);
+               break;
+         }
+         break;
+      }
+      case 1:
+      case 2:
+      {
+         // Clockwise rotation in 2D around the origin
+         const double w = M_PI/2;
+         switch (dim)
+         {
+            case 1: v(0) = 1.0; break;
+            case 2: v(0) = w*X(1); v(1) = -w*X(0); break;
+            case 3: v(0) = w*X(1); v(1) = -w*X(0); v(2) = 0.0; break;
+         }
+         break;
+      }
+      case 3:
+      {
+         // Clockwise twisting rotation in 2D around the origin
+         const double w = M_PI/2;
+         double d = max((X(0)+1.)*(1.-X(0)),0.) * max((X(1)+1.)*(1.-X(1)),0.);
+         d = d*d;
+         switch (dim)
+         {
+            case 1: v(0) = 1.0; break;
+            case 2: v(0) = d*w*X(1); v(1) = -d*w*X(0); break;
+            case 3: v(0) = d*w*X(1); v(1) = -d*w*X(0); v(2) = 0.0; break;
+         }
+         break;
+      }
+   }
+}
+
+// Initial condition
+double u0_function(const Vector &x)
+{
+   int dim = x.Size();
+
+   // map to the reference [-1,1] domain
+   Vector X(dim);
+   for (int i = 0; i < dim; i++)
+   {
+      double center = (bb_min[i] + bb_max[i]) * 0.5;
+      X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
+   }
+
+   switch (problem)
+   {
+      case 0:
+      case 1:
+      {
+         switch (dim)
+         {
+            case 1:
+               return exp(-40.*pow(X(0)-0.5,2));
+            case 2:
+            case 3:
+            {
+               double rx = 0.45, ry = 0.25, cx = 0., cy = -0.2, w = 10.;
+               if (dim == 3)
+               {
+                  const double s = (1. + 0.25*cos(2*M_PI*X(2)));
+                  rx *= s;
+                  ry *= s;
+               }
+               return ( erfc(w*(X(0)-cx-rx))*erfc(-w*(X(0)-cx+rx)) *
+                        erfc(w*(X(1)-cy-ry))*erfc(-w*(X(1)-cy+ry)) )/16 + 1e-09;
+            }
+         }
+      }
+      case 2:
+      {
+         double x_ = X(0), y_ = X(1), rho, phi;
+         rho = hypot(x_, y_);
+         phi = atan2(y_, x_);
+         return pow(sin(M_PI*rho),2)*sin(3*phi);
+      }
+      case 3:
+      {
+         const double f = M_PI;
+         return sin(f*X(0))*sin(f*X(1));
+      }
+   }
+   return 0.0;
+}
+
+// Inflow boundary condition (zero for the problems considered in this example)
+double inflow_function(const Vector &x)
+{
+   switch (problem)
+   {
+      case 0:
+      case 1:
+      case 2:
+      case 3: return 0.0;
+   }
+   return 0.0;
+}

examples/proxGalerkinAdvBlockOp.cpp

@@ -0,0 +1,794 @@
+//                                MFEM Example 9
+//
+// Compile with: make ex9
+//
+// Sample runs:
+//    ex9 -m ../data/periodic-segment.mesh -p 0 -r 2 -dt 0.005
+//    ex9 -m ../data/periodic-square.mesh -p 0 -r 2 -dt 0.01 -tf 10
+//    ex9 -m ../data/periodic-hexagon.mesh -p 0 -r 2 -dt 0.01 -tf 10
+//    ex9 -m ../data/periodic-square.mesh -p 1 -r 2 -dt 0.005 -tf 9
+//    ex9 -m ../data/periodic-hexagon.mesh -p 1 -r 2 -dt 0.005 -tf 9
+//    ex9 -m ../data/amr-quad.mesh -p 1 -r 2 -dt 0.002 -tf 9
+//    ex9 -m ../data/amr-quad.mesh -p 1 -r 2 -dt 0.02 -s 13 -tf 9
+//    ex9 -m ../data/star-q3.mesh -p 1 -r 2 -dt 0.005 -tf 9
+//    ex9 -m ../data/star-mixed.mesh -p 1 -r 2 -dt 0.005 -tf 9
+//    ex9 -m ../data/disc-nurbs.mesh -p 1 -r 3 -dt 0.005 -tf 9
+//    ex9 -m ../data/disc-nurbs.mesh -p 2 -r 3 -dt 0.005 -tf 9
+//    ex9 -m ../data/periodic-square.mesh -p 3 -r 4 -dt 0.0025 -tf 9 -vs 20
+//    ex9 -m ../data/periodic-cube.mesh -p 0 -r 2 -o 2 -dt 0.02 -tf 8
+//    ex9 -m ../data/periodic-square.msh -p 0 -r 2 -dt 0.005 -tf 2
+//    ex9 -m ../data/periodic-cube.msh -p 0 -r 1 -o 2 -tf 2
+//
+// Device sample runs:
+//    ex9 -pa
+//    ex9 -ea
+//    ex9 -fa
+//    ex9 -pa -m ../data/periodic-cube.mesh
+//    ex9 -pa -m ../data/periodic-cube.mesh -d cuda
+//    ex9 -ea -m ../data/periodic-cube.mesh -d cuda
+//    ex9 -fa -m ../data/periodic-cube.mesh -d cuda
+//    ex9 -pa -m ../data/amr-quad.mesh -p 1 -r 2 -dt 0.002 -tf 9 -d cuda
+//
+// Description:  This example code solves the time-dependent advection equation
+//               du/dt + v.grad(u) = 0, where v is a given fluid velocity, and
+//               u0(x)=u(0,x) is a given initial condition.
+//
+//               The example demonstrates the use of Discontinuous Galerkin (DG)
+//               bilinear forms in MFEM (face integrators), the use of implicit
+//               and explicit ODE time integrators, the definition of periodic
+//               boundary conditions through periodic meshes, as well as the use
+//               of GLVis for persistent visualization of a time-evolving
+//               solution. The saving of time-dependent data files for external
+//               visualization with VisIt (visit.llnl.gov) and ParaView
+//               (paraview.org) is also illustrated.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+#include <algorithm>
+
+using namespace std;
+using namespace mfem;
+
+// Choice for the problem setup. The fluid velocity, initial condition and
+// inflow boundary condition are chosen based on this parameter.
+
+// Velocity coefficient
+std::function<void(const Vector&, Vector&)> getVelocityFunction(int problem);
+// Initial condition
+std::function<double(const Vector&)> getInflowFunction(int problem);
+
+// Inflow boundary condition
+std::function<double(const Vector&)> getInitialCondition(int problem);
+
+// Mesh bounding box
+Vector bb_min, bb_max;
+
+#ifdef MFEM_USE_LAPACK
+/// Integrator that inverts the matrix assembled by another integrator via Pseudo-inverse.
+class PseudoInverseIntegrator : public BilinearFormIntegrator
+{
+private:
+   int own_integrator;
+   BilinearFormIntegrator *integrator;
+
+public:
+   PseudoInverseIntegrator(BilinearFormIntegrator *integ, int own_integ = 1)
+   { integrator = integ; own_integrator = own_integ; }
+
+   virtual void SetIntRule(const IntegrationRule *ir);
+
+   virtual void AssembleElementMatrix(const FiniteElement &el,
+                                      ElementTransformation &Trans,
+                                      DenseMatrix &elmat);
+
+   virtual ~PseudoInverseIntegrator() { if (own_integrator) { delete integrator; } }
+};
+#endif
+
+class DG_Solver : public Solver
+{
+private:
+   SparseMatrix &M, &K, A;
+   GMRESSolver linear_solver;
+   BlockILU prec;
+   double dt;
+public:
+   DG_Solver(SparseMatrix &M_, SparseMatrix &K_, const FiniteElementSpace &fes)
+      : M(M_),
+        K(K_),
+        prec(fes.GetFE(0)->GetDof(),
+             BlockILU::Reordering::MINIMUM_DISCARDED_FILL),
+        dt(-1.0)
+   {
+      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);
+   }
+
+   void SetTimeStep(double dt_)
+   {
+      if (dt_ != dt)
+      {
+         dt = dt_;
+         // Form operator A = M - dt*K
+         A = K;
+         A *= -dt;
+         A += M;
+
+         // this will also call SetOperator on the preconditioner
+         linear_solver.SetOperator(A);
+      }
+   }
+
+   void SetOperator(const Operator &op)
+   {
+      linear_solver.SetOperator(op);
+   }
+
+   virtual void Mult(const Vector &x, Vector &y) const
+   {
+      linear_solver.Mult(x, y);
+   }
+};
+
+/** 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
+    be written as a general ODE, du/dt = M^{-1} (K u + b), and this class is
+    used to evaluate the right-hand side. */
+class FE_Evolution : public TimeDependentOperator
+{
+private:
+   BilinearForm &M, &K;
+   const Vector &b;
+   Solver *M_prec;
+   CGSolver M_solver;
+   DG_Solver *dg_solver;
+   FiniteElementSpace *fes;
+
+   Array<int> offsets;
+   BlockVector sol;
+   Vector ktemp;
+   GridFunction latent;
+   GridFunction latent_k;
+   GridFunction xnew;
+
+   GridFunctionCoefficient latent_cf;
+   GridFunctionCoefficient latent_k_cf;
+   GridFunctionCoefficient xnew_cf;
+   Vector delta_k, delta_latent, delta_latent_inner;
+
+   TransformedCoefficient expLatent_cf;
+   SumCoefficient diff_latent_cf;
+   TransformedCoefficient expResidual_cf;
+
+   LinearForm delta_latent_form;
+   NonlinearForm invExpLatentForm;
+   LinearForm newtonRHS;
+
+
+   mutable Vector z;
+
+public:
+   FE_Evolution(BilinearForm &M_, BilinearForm &K_, const Vector &b_);
+
+   virtual void Mult(const Vector &x, Vector &y) const;
+   virtual void ImplicitSolve(const double dt, const Vector &x, Vector &k);
+
+   virtual ~FE_Evolution();
+   bool postprocess = false;
+   bool initialized;
+};
+
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   int problem = 1;
+   const char *mesh_file = "../data/periodic-square-4x4.mesh";
+   int ref_levels = 2;
+   int order = 3;
+   bool pa = false;
+   bool ea = false;
+   bool fa = false;
+   const char *device_config = "cpu";
+   int ode_solver_type = 11;
+   double t_final = 2;
+   double dt = 0.01;
+   bool visualization = true;
+   bool visit = false;
+   bool paraview = false;
+   bool binary = false;
+   int vis_steps = 1;
+   bool applyPostprocess = false;
+
+   int precision = 8;
+   cout.precision(precision);
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&problem, "-p", "--problem",
+                  "Problem setup to use. See options in velocity_function().");
+   args.AddOption(&ref_levels, "-r", "--refine",
+                  "Number of times to refine the mesh uniformly.");
+   args.AddOption(&order, "-o", "--order",
+                  "Order (degree) of the finite elements.");
+   args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
+                  "--no-partial-assembly", "Enable Partial Assembly.");
+   args.AddOption(&ea, "-ea", "--element-assembly", "-no-ea",
+                  "--no-element-assembly", "Enable Element Assembly.");
+   args.AddOption(&fa, "-fa", "--full-assembly", "-no-fa",
+                  "--no-full-assembly", "Enable Full Assembly.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&ode_solver_type, "-s", "--ode-solver",
+                  "ODE solver: 1 - Forward Euler,\n\t"
+                  "            2 - RK2 SSP, 3 - RK3 SSP, 4 - RK4, 6 - RK6,\n\t"
+                  "            11 - Backward Euler,\n\t"
+                  "            12 - SDIRK23 (L-stable), 13 - SDIRK33,\n\t"
+                  "            22 - Implicit Midpoint Method,\n\t"
+                  "            23 - SDIRK23 (A-stable), 24 - SDIRK34");
+   args.AddOption(&t_final, "-tf", "--t-final",
+                  "Final time; start time is 0.");
+   args.AddOption(&dt, "-dt", "--time-step",
+                  "Time step.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&visit, "-visit", "--visit-datafiles", "-no-visit",
+                  "--no-visit-datafiles",
+                  "Save data files for VisIt (visit.llnl.gov) visualization.");
+   args.AddOption(&paraview, "-paraview", "--paraview-datafiles", "-no-paraview",
+                  "--no-paraview-datafiles",
+                  "Save data files for ParaView (paraview.org) visualization.");
+   args.AddOption(&binary, "-binary", "--binary-datafiles", "-ascii",
+                  "--ascii-datafiles",
+                  "Use binary (Sidre) or ascii format for VisIt data files.");
+   args.AddOption(&vis_steps, "-vs", "--visualization-steps",
+                  "Visualize every n-th timestep.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   Device device(device_config);
+   device.Print();
+
+   // 2. Read the mesh from the given mesh file. We can handle geometrically
+   //    periodic meshes in this code.
+   Mesh mesh(mesh_file, 1, 1);
+   int dim = mesh.Dimension();
+
+   // 3. Define the ODE solver used for time integration. Several explicit
+   //    Runge-Kutta methods are available.
+   ODESolver *ode_solver = NULL;
+   switch (ode_solver_type)
+   {
+      // Explicit methods
+      case 1: ode_solver = new ForwardEulerSolver; break;
+      case 2: ode_solver = new RK2Solver(1.0); break;
+      case 3: ode_solver = new RK3SSPSolver; break;
+      case 4: ode_solver = new RK4Solver; break;
+      case 6: ode_solver = new RK6Solver; break;
+      // Implicit (L-stable) methods
+      case 11: ode_solver = new BackwardEulerSolver; break;
+      case 12: ode_solver = new SDIRK23Solver(2); break;
+      case 13: ode_solver = new SDIRK33Solver; break;
+      // Implicit A-stable methods (not L-stable)
+      case 22: ode_solver = new ImplicitMidpointSolver; break;
+      case 23: ode_solver = new SDIRK23Solver; break;
+      case 24: ode_solver = new SDIRK34Solver; break;
+
+      default:
+         cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
+         return 3;
+   }
+
+   // 4. Refine the mesh to increase the resolution. In this example we do
+   //    'ref_levels' of uniform refinement, where 'ref_levels' is a
+   //    command-line parameter. If the mesh is of NURBS type, we convert it to
+   //    a (piecewise-polynomial) high-order mesh.
+   for (int lev = 0; lev < ref_levels; lev++)
+   {
+      mesh.UniformRefinement();
+   }
+   if (mesh.NURBSext)
+   {
+      mesh.SetCurvature(max(order, 1));
+   }
+   mesh.GetBoundingBox(bb_min, bb_max, max(order, 1));
+
+   // 5. Define the discontinuous DG finite element space of the given
+   //    polynomial order on the refined mesh.
+   DG_FECollection fec(order, dim, BasisType::GaussLobatto);
+   FiniteElementSpace fes(&mesh, &fec);
+
+   cout << "Number of unknowns: " << fes.GetVSize() << endl;
+
+   // 6. Set up and assemble the bilinear and linear forms corresponding to the
+   //    DG discretization. The DGTraceIntegrator involves integrals over mesh
+   //    interior faces.
+   VectorFunctionCoefficient velocity(dim, getVelocityFunction(problem));
+   FunctionCoefficient inflow(getInflowFunction(problem));
+   FunctionCoefficient u0(getInitialCondition(problem));
+
+   BilinearForm m(&fes);
+   BilinearForm k(&fes);
+   if (pa)
+   {
+      m.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+      k.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+   }
+   else if (ea)
+   {
+      m.SetAssemblyLevel(AssemblyLevel::ELEMENT);
+      k.SetAssemblyLevel(AssemblyLevel::ELEMENT);
+   }
+   else if (fa)
+   {
+      m.SetAssemblyLevel(AssemblyLevel::FULL);
+      k.SetAssemblyLevel(AssemblyLevel::FULL);
+   }
+   m.AddDomainIntegrator(new MassIntegrator);
+   constexpr double alpha = -1.0;
+   k.AddDomainIntegrator(new ConvectionIntegrator(velocity, alpha));
+   k.AddInteriorFaceIntegrator(
+      new NonconservativeDGTraceIntegrator(velocity, alpha));
+   k.AddBdrFaceIntegrator(
+      new NonconservativeDGTraceIntegrator(velocity, alpha));
+
+   LinearForm b(&fes);
+   b.AddBdrFaceIntegrator(
+      new BoundaryFlowIntegrator(inflow, velocity, alpha));
+
+   m.Assemble();
+   int skip_zeros = 0;
+   k.Assemble(skip_zeros);
+   b.Assemble();
+   m.Finalize();
+   k.Finalize(skip_zeros);
+
+   // 7. Define the initial conditions, save the corresponding grid function to
+   //    a file and (optionally) save data in the VisIt format and initialize
+   //    GLVis visualization.
+   GridFunction u(&fes);
+   u.ProjectCoefficient(u0);
+
+   {
+      ofstream omesh("ex9.mesh");
+      omesh.precision(precision);
+      mesh.Print(omesh);
+      ofstream osol("ex9-init.gf");
+      osol.precision(precision);
+      u.Save(osol);
+   }
+
+   // Create data collection for solution output: either VisItDataCollection for
+   // ascii data files, or SidreDataCollection for binary data files.
+   DataCollection *dc = NULL;
+   if (visit)
+   {
+      if (binary)
+      {
+#ifdef MFEM_USE_SIDRE
+         dc = new SidreDataCollection("Example9", &mesh);
+#else
+         MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
+#endif
+      }
+      else
+      {
+         dc = new VisItDataCollection("Example9", &mesh);
+         dc->SetPrecision(precision);
+      }
+      dc->RegisterField("solution", &u);
+      dc->SetCycle(0);
+      dc->SetTime(0.0);
+      dc->Save();
+   }
+
+   ParaViewDataCollection *pd = NULL;
+   if (paraview)
+   {
+      pd = new ParaViewDataCollection("Example9", &mesh);
+      pd->SetPrefixPath("ParaView");
+      pd->RegisterField("solution", &u);
+      pd->SetLevelsOfDetail(order);
+      pd->SetDataFormat(VTKFormat::BINARY);
+      pd->SetHighOrderOutput(true);
+      pd->SetCycle(0);
+      pd->SetTime(0.0);
+      pd->Save();
+   }
+
+   socketstream sout;
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      sout.open(vishost, visport);
+      if (!sout)
+      {
+         cout << "Unable to connect to GLVis server at "
+              << vishost << ':' << visport << endl;
+         visualization = false;
+         cout << "GLVis visualization disabled.\n";
+      }
+      else
+      {
+         sout.precision(precision);
+         sout << "solution\n" << mesh << u;
+         // sout << "pause\n";
+         sout << flush;
+         // cout << "GLVis visualization paused."
+         //      << " Press space (in the GLVis window) to resume it.\n";
+      }
+   }
+
+   // 8. 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);
+   adv.postprocess = applyPostprocess;
+
+   double t = 0.0;
+   adv.SetTime(t);
+   ode_solver->Init(adv);
+
+   bool done = false;
+   for (int ti = 0; !done; )
+   {
+      double dt_real = min(dt, t_final - t);
+      ode_solver->Step(u, t, dt_real);
+      ti++;
+
+      done = (t >= t_final - 1e-8*dt);
+
+      if (done || ti % vis_steps == 0)
+      {
+         cout << "time step: " << ti << ", time: " << t << endl;
+
+         if (visualization)
+         {
+            sout << "solution\n" << mesh << u << flush;
+         }
+
+         if (visit)
+         {
+            dc->SetCycle(ti);
+            dc->SetTime(t);
+            dc->Save();
+         }
+
+         if (paraview)
+         {
+            pd->SetCycle(ti);
+            pd->SetTime(t);
+            pd->Save();
+         }
+      }
+   }
+
+   // 9. Save the final solution. This output can be viewed later using GLVis:
+   //    "glvis -m ex9.mesh -g ex9-final.gf".
+   {
+      ofstream osol("ex9-final.gf");
+      osol.precision(precision);
+      u.Save(osol);
+   }
+   out << "L2 error: " << u.ComputeL2Error(u0) << std::endl;
+
+   // 10. Free the used memory.
+   delete ode_solver;
+   delete pd;
+   delete dc;
+
+   return 0;
+}
+
+
+// Implementation of class FE_Evolution
+FE_Evolution::FE_Evolution(BilinearForm &M_, BilinearForm &K_, const Vector &b_)
+   : TimeDependentOperator(M_.Height()), M(M_), K(K_), b(b_), fes(M.FESpace()),
+     offsets({0, fes->GetNDofs(), fes->GetNDofs()*2}), sol(offsets),
+ktemp(sol.GetBlock(0).GetData(), fes->GetNDofs()),
+latent(fes, sol.GetBlock(1).GetData()), latent_k(fes), xnew(fes),
+latent_cf(&latent),
+latent_k_cf(&latent_k), xnew_cf(&xnew), delta_k(xnew.Size()),
+delta_latent(latent.Size()), delta_latent_inner(latent.Size()),
+expLatent_cf(&latent_cf, [](double x) {return exp(x);}),
+diff_latent_cf(latent_cf, latent_k_cf, -1.0, 1.0),
+expResidual_cf(&latent_cf, [](double psi) {return exp(psi)*(1.0 - psi);}),
+delta_latent_form(fes), invExpLatentForm(fes), newtonRHS(fes), z(M_.Height()),
+initialized(false)
+{
+
+   delta_latent_form.AddDomainIntegrator(new DomainLFIntegrator(diff_latent_cf));
+   // since DG is local space, we can perform inverse element-wise.
+   invExpLatentForm.AddDomainIntegrator(new InverseIntegrator(new MassIntegrator(
+                                                                 expLatent_cf)));
+   newtonRHS.AddDomainIntegrator(new DomainLFIntegrator(expResidual_cf));
+   Array<int> ess_tdof_list;
+   if (M.GetAssemblyLevel() == AssemblyLevel::LEGACY)
+   {
+      M_prec = new DSmoother(M.SpMat());
+      M_solver.SetOperator(M.SpMat());
+      dg_solver = new DG_Solver(M.SpMat(), K.SpMat(), *M.FESpace());
+   }
+   else
+   {
+      M_prec = new OperatorJacobiSmoother(M, ess_tdof_list);
+      M_solver.SetOperator(M);
+      dg_solver = NULL;
+   }
+   M_solver.SetPreconditioner(*M_prec);
+   M_solver.iterative_mode = false;
+   M_solver.SetRelTol(1e-9);
+   M_solver.SetAbsTol(0.0);
+   M_solver.SetMaxIter(100);
+   M_solver.SetPrintLevel(0);
+}
+
+void FE_Evolution::Mult(const Vector &x, Vector &y) const
+{
+   // y = M^{-1} (K x + b)
+   K.Mult(x, z);
+   z += b;
+   M_solver.Mult(z, y);
+}
+
+void FE_Evolution::ImplicitSolve(const double dt, const Vector &x, Vector &k)
+{
+   if (!initialized)
+   {
+      ktemp = 0.0;
+      latent = x;
+      // latent.ApplyMap([](double x) {return log(x); });
+      latent.ApplyMap([](double x) {return x > 1e-09 ? log(x) : -20.0; });
+      initialized = true;
+   }
+   int printlevel = 2;
+   bool debug = true;
+   MFEM_VERIFY(dg_solver != NULL,
+               "Implicit time integration is not supported with partial assembly");
+   double alpha_prev = 1.0;
+   // auto alphamaker = [&alpha_prev](int i) {return min(max(1.0, pow(1.5, pow(1.5, i-1)) - alpha_prev), 1e10);};
+   auto alphamaker = [](int i) {return (double)(i*i);};
+   GridFunction zero_gf(fes);
+   zero_gf = 0.0;
+
+   BlockOperator op(offsets);
+   SparseMatrix A(x.Size());
+   A = K.SpMat();
+   A *= -dt;
+   A += M.SpMat();
+   op.SetBlock(0, 0, &A);
+
+   SparseMatrix offdiag(x.Size());
+   offdiag = M.SpMat();
+   offdiag *= dt;
+
+   op.SetBlock(1, 0, &offdiag);
+   op.SetBlock(0, 1, &offdiag);
+
+   BlockVector rhs(offsets);
+   z.SetDataAndSize(rhs.GetBlock(0).GetData(), x.Size());
+   Vector z2(rhs.GetBlock(1).GetData(), x.Size());
+
+
+   BlockDiagonalPreconditioner prec(offsets);
+   BilinearForm invM(fes);
+   invM.AddDomainIntegrator(new InverseIntegrator(new MassIntegrator()));
+   invM.Assemble();
+   prec.SetDiagonalBlock(0, &invM.SpMat());
+
+   int maxIter(100000);
+   double rtol(1.e-6);
+   double atol(1.e-10);
+   GMRESSolver solver;
+   solver.SetAbsTol(atol);
+   solver.SetRelTol(rtol);
+   solver.SetMaxIter(maxIter);
+   solver.SetOperator(op);
+   solver.SetPreconditioner(prec);
+   solver.SetPrintLevel(0);
+   solver.iterative_mode = true;
+   bool converged = false;
+   for (int i=1; i<1e04; i++)
+   {
+      latent_k = latent;
+      if (i > 1)
+      {
+         double alpha = alphamaker(i);
+         double dalpha = alpha / alpha_prev;
+         alpha_prev = alpha;
+
+         A *= dalpha;
+         invM.SpMat().operator*=(dalpha);
+      }
+      bool converged_inner = false;
+      for (int j=1; j<=10; j++)
+      {
+         delta_k = ktemp;
+         rhs = 0.0;
+
+         K.Mult(x, z);
+         z += b;
+         z *= alphamaker(i);
+         M.AddMult(latent, z, -dt);
+         M.AddMult(latent_k, z, dt);
+
+         BilinearForm expM(fes);
+         expM.AddDomainIntegrator(new MassIntegrator(expLatent_cf));
+         expM.Assemble();
+         expM.SpMat().operator*=(-1.0);
+
+         op.SetBlock(1, 1, &expM.SpMat());
+
+         if (debug)
+         {
+            out << latent.Min() << ", " << latent.Max() << std::endl;
+            if (latent.CheckFinite())
+            {
+               for (auto val : latent) {if (!IsFinite(val)) out << val << " .. ";}
+               out << std::endl;
+               mfem_error("latent variable is not finite.");
+            }
+            GridFunction dummy(fes);
+            dummy.ProjectCoefficient(expLatent_cf);
+            if (dummy.CheckFinite())
+            {
+               for (auto val : dummy) {if (!IsFinite(val)) out << val << " .. ";}
+               out << std::endl;
+               mfem_error("Mapped latent variable is not finite.");
+            }
+         }
+         BilinearForm expInvM(fes);
+
+         expInvM.AddDomainIntegrator(new PseudoInverseIntegrator(new MassIntegrator(
+                                                                    expLatent_cf)));
+         expInvM.Assemble();
+         expInvM.SpMat().operator*=(-1.0);
+         prec.SetDiagonalBlock(1, &expInvM.SpMat());
+
+         M.Mult(x, z2);
+         z2.Neg();
+         newtonRHS.Assemble();
+         z2 += newtonRHS;
+
+         solver.Mult(rhs, sol);
+         if (!solver.GetConverged())
+         {
+            mfem_warning("LinearSolver failed to converge.");
+         }
+
+         delta_k -= ktemp;
+         if (printlevel > 1) { out << "\t\tInner Step " << j << ", " << delta_k.Normlinf() << std::endl; }
+         if (delta_k.Normlinf() < 1e-07)
+         {
+            if (printlevel > 0) { out << "\t\tInnerloop convereged in " << j << " steps."<< std::endl; }
+            converged_inner = true;
+            break;
+         }
+      }
+      if (!converged_inner) {mfem_warning("Subproblem failed to converge."); }
+      TransformedCoefficient latent_diff(&latent_k_cf, &latent_cf, [](double x,
+      double xnew) {return exp(x)*(x-xnew);});
+      double residual = zero_gf.ComputeLpError(infinity(), latent_diff);
+      if (printlevel > 1) { out << "\tOuter Step " << i << ", " << residual << std::endl; }
+      if (zero_gf.ComputeLpError(infinity(), latent_diff) < 1e-07)
+      {
+         converged = true;
+         if (printlevel > 0) out << "\tOuterloop convereged in " << i <<
+                                    " steps.\n---------------------------------------------------"<< std::endl;
+         break;
+      }
+   }
+   if (!converged) {mfem_warning("Main loop failed to converge."); }
+   if (postprocess)
+   {
+      k = latent;
+      k.ApplyMap([](double x) {return exp(x); });
+      k -= x;
+      k *= dt;
+   }
+   else
+   {
+      k = ktemp;
+   }
+}
+
+FE_Evolution::~FE_Evolution()
+{
+   delete M_prec;
+   delete dg_solver;
+}
+
+// Velocity coefficient
+std::function<void(const Vector&, Vector &)> getVelocityFunction(int problem)
+{
+   switch (problem)
+   {
+      case 0:
+         return [](const Vector &x, Vector &y) { y = 1.0; };
+      case 1:
+         return [](const Vector &x, Vector &y) { y = 1.0; };
+      default:
+         mfem_error("Undefined problem");
+         return [](const Vector &x, Vector &y) { };
+   }
+
+}
+
+// Initial condition
+std::function<double(const Vector &)> getInitialCondition(int problem)
+{
+   switch (problem)
+   {
+      case 0:
+         return [](const Vector &x)
+         {
+            double rx = 0.45, ry = 0.25, cx = 0., cy = -0.2, w = 10.;
+            return ( erfc(w*(x(0)-cx-rx))*erfc(-w*(x(0)-cx+rx)) *
+                     erfc(w*(x(1)-cy-ry))*erfc(-w*(x(1)-cy+ry)) )/16;
+         };
+      case 1:
+         return [](const Vector &x)
+         {
+            double offset = 1e-04; // away from zero
+            double y = 1.0;
+            for (auto v:x) {y *= sin(M_PI*v); }
+            y += 1.0 + offset;
+            return y;
+         };
+      default:
+         mfem_error("Undefined problem");
+         return [](const Vector &x) { return 0.0; };
+   }
+}
+
+// Inflow boundary condition (zero for the problems considered in this example)
+std::function<double(const Vector&)> getInflowFunction(int problem)
+{
+   switch (problem)
+   {
+      case 0:
+         return [](const Vector &x) {return 0.0; };
+      case 1:
+         return [](const Vector &x) {return 0.0; };
+      default:
+         mfem_error("Undefined problem");
+         return [](const Vector &x) { return 0.0; };
+   }
+}
+
+
+#ifdef MFEM_USE_LAPACK
+void PseudoInverseIntegrator::SetIntRule(const IntegrationRule *ir)
+{
+   IntRule = ir;
+   integrator->SetIntRule(ir);
+}
+
+void PseudoInverseIntegrator::AssembleElementMatrix(
+   const FiniteElement &el, ElementTransformation &Trans, DenseMatrix &elmat)
+{
+   integrator->AssembleElementMatrix(el, Trans, elmat);
+   DenseMatrixSVD svd(elmat,true,true);
+   svd.Eval(elmat);
+   Vector &sigma = svd.Singularvalues();
+   DenseMatrix &U = svd.LeftSingularvectors();
+   DenseMatrix &V = svd.RightSingularvectors();
+   sigma.ApplyMap([](double x) {return x > 1e-08 ? 1.0 / x : 0; });
+   DenseMatrix Vt(V); Vt.Transpose();
+   MultADBt(U,sigma,Vt,elmat);
+}
+#endif

examples/test_RT.cpp

@@ -0,0 +1,33 @@
+#include "mfem.hpp"
+#include <chrono>
+
+
+using namespace mfem;
+int main()
+{
+   const int order = 2;
+   const int dim = 2;
+   for (int i=0; i< 5; i++)
+   {
+      Mesh mesh = Mesh::MakeCartesian2D(static_cast<int>(std::pow(2, i)),
+                                        static_cast<int>(std::pow(2, i)), Element::Type::QUADRILATERAL);
+      mesh.EnsureNCMesh();
+      out<<mesh.GetNE() <<std::endl;
+
+      FiniteElementCollection *fec = new RT_FECollection(order, dim);
+      FiniteElementSpace *fes = new FiniteElementSpace(&mesh, fec);
+
+      VectorFunctionCoefficient v(2, [](const Vector &x, double t, Vector &y)
+      {
+         y(0) = 1.0;
+         y(1) = 0.0;
+      });
+
+      GridFunction x(fes);
+      x.ProjectCoefficient(v);
+      out << "(" << x.Min() << ", " << x.Max() << ")" << std::endl;
+      delete fes;
+      delete fec;
+   }
+   return 0;
+}

examples/test_estimator.cpp

@@ -0,0 +1,61 @@
+#include "mfem.hpp"
+#include <chrono>
+
+
+using namespace mfem;
+int main()
+{
+   const int order = 2;
+   const int dim = 2;
+   for (int i=0; i<12; i++)
+   {
+      Mesh mesh = Mesh::MakeCartesian2D(static_cast<int>(std::pow(2, i)),
+                                        static_cast<int>(std::pow(2, i)), Element::Type::QUADRILATERAL);
+      mesh.EnsureNCMesh();
+      out<<mesh.GetNE() <<std::endl;
+
+      FiniteElementCollection *fec = new DG_FECollection(order, dim);
+      FiniteElementSpace *fes = new FiniteElementSpace(&mesh, fec, 2,
+                                                       Ordering::byNODES);
+
+      VectorFunctionCoefficient v(2, [](const Vector &x, double t, Vector &y)
+      {
+         y(0) = std::sin(x(0))*std::sin(x(1));
+         y(1) = std::cos(x(0))*std::cos(x(1));
+      });
+
+      GridFunction x(fes);
+      x.ProjectCoefficient(v);
+
+      out << "Estimation" << std::endl;
+      auto start_1 = std::chrono::high_resolution_clock::now();
+      ProjectionErrorEstimator estimator(x);
+      auto & estimators = estimator.GetLocalErrors();
+      out << "estimated error: " << estimator.GetTotalError() << std::endl;
+      auto stop_1 = std::chrono::high_resolution_clock::now();
+      auto duration_1 = std::chrono::duration_cast<std::chrono::microseconds>
+                        (stop_1 - start_1);
+      out << "wall time: " << duration_1.count() / 1000.0 << "ms" << std::endl;
+
+      out << "Estimation" << std::endl;
+      auto start_2 = std::chrono::high_resolution_clock::now();
+      PRefDiffEstimator estimator_Katen(x, -1);
+      auto & estimators_Katen = estimator_Katen.GetLocalErrors();
+      out << "estimated error: " << estimator_Katen.GetTotalError() << std::endl;
+      auto stop_2 = std::chrono::high_resolution_clock::now();
+      auto duration_2 = std::chrono::duration_cast<std::chrono::microseconds>
+                        (stop_2 - start_2);
+      out << "wall time: " << duration_2.count() / 1000.0 <<  "ms" << std::endl;
+
+      // double total_error = estimators.Norml2();
+      out << x.ComputeL2Error(v) << std::endl;
+      out << estimator.GetTotalError() / x.ComputeL2Error(v) << std::endl;
+      out << estimator_Katen.GetTotalError() / x.ComputeL2Error(v) << std::endl;
+      // out << estimators.Max() << std::endl;
+      delete fes;
+      delete fec;
+   }
+
+   return 0;
+
+}

fem/CMakeLists.txt

@@ -146,6 +146,7 @@ set(SRCS
   tmop_amr.cpp
   gslib.cpp
   transfer.cpp
+  hyperbolic_conservation_laws.cpp
   )
 
 set(HDRS
@@ -237,6 +238,7 @@ set(HDRS
   tmop_amr.hpp
   gslib.hpp
   transfer.hpp
+  hyperbolic_conservation_laws.hpp
   )
 
 if (MFEM_USE_SIDRE)

fem/estimators.cpp

@@ -11,6 +11,8 @@
 
 #include "estimators.hpp"
 
+#include "transfer.hpp"
+
 namespace mfem
 {
 
@@ -21,25 +23,120 @@ void ZienkiewiczZhuEstimator::ComputeEstimates()
    // ParFiniteElementSpace and 'solution' is a ParGridFunction.
    GridFunction flux(flux_space);
 
-   if (!anisotropic) { aniso_flags.SetSize(0); }
+   if (!anisotropic)
+   {
+      aniso_flags.SetSize(0);
+   }
    total_error = ZZErrorEstimator(integ, solution, flux, error_estimates,
                                   anisotropic ? &aniso_flags : NULL,
-                                  flux_averaging,
-                                  with_coeff);
+                                  flux_averaging, with_coeff);
 
    current_sequence = solution.FESpace()->GetMesh()->GetSequence();
 }
 
+PRefDiffEstimator::PRefDiffEstimator(GridFunction &sol_, int p_comp_)
+   : solution(&sol_), p_comp(p_comp_) {}
+
+void PRefDiffEstimator::ComputeEstimates()
+{
+   const int nelem = solution->FESpace()->GetNE();
+   error_estimates.SetSize(nelem);
+
+   FiniteElementSpace *fespace = solution->FESpace();
+   FiniteElementSpace fespaceComp = FiniteElementSpace(*(solution->FESpace()));
+
+   for (int e = 0; e < nelem; e++)
+   {
+      fespaceComp.SetElementOrder(e, solution->FESpace()->GetElementOrder(e));
+   }
+   fespaceComp.Update(false);
+
+   GridFunction solutionComp(&fespaceComp);
+   solutionComp = *solution;
+
+   for (int e = 0; e < nelem; e++)
+   {
+      int setOrder = p_comp >= 0
+                     ? p_comp
+                     : std::max(0, fespace->GetElementOrder(e) + p_comp);
+      fespaceComp.SetElementOrder(e, setOrder);
+   }
+   fespaceComp.Update(false);
+   solutionComp.Update();
+
+   PRefinementTransferOperator Transfer(*fespace, fespaceComp);
+   Transfer.Mult(*solution, solutionComp);
+
+   GridFunctionCoefficient solutionCompCoeff(&solutionComp);
+   solution->ComputeElementL2Errors(solutionCompCoeff, error_estimates);
+
+   total_error = error_estimates.Norml2();
+}
+
+void ProjectionErrorEstimator::ComputeEstimates()
+{
+   // Gather information and make projection space
+   const int nelem = solution.FESpace()->GetNE();     // number of elements
+   FiniteElementSpace *fespace = solution.FESpace();  // solution fespace
+   FiniteElementSpace projectionSpace(*fespace);
+   int max_order = 0, min_order = INT_MAX;
+   for (int i = 0; i < nelem; i++)
+   {
+      int order = std::max(0, fespace->GetElementOrder(i) - offset);
+      projectionSpace.SetElementOrder(i, order);  // update polynomial order
+   }
+   projectionSpace.Update(false);
+
+   // local mass inverse
+   DenseMatrix invMe;  // auxiliary local mass matrix
+   InverseIntegrator invmi(new MassIntegrator());
+
+   // (u_h, v)
+   VectorGridFunctionCoefficient sol_gf(&solution);
+   VectorDomainLFIntegrator int_sol(sol_gf);
+
+   // Projected Space
+   GridFunction projectedSol(&projectionSpace);
+   Array<int> dofs;
+   Vector loc_proj, rhs;
+
+   // resize it to the current number of elements
+   for (int i = 0; i < projectionSpace.GetNE(); i++)
+   {
+      int nDofs = projectionSpace.GetFE(i)->GetDof();
+      dofs.SetSize(nDofs);
+      invMe.SetSize(nDofs);
+      loc_proj.SetSize(nDofs * projectionSpace.GetVDim());
+
+      projectionSpace.GetElementVDofs(i, dofs);
+      invmi.AssembleElementMatrix(*projectionSpace.GetFE(i),
+                                  *projectionSpace.GetElementTransformation(i),
+                                  invMe);
+      rhs.SetSize(nDofs * projectionSpace.GetVDim());
+      int_sol.AssembleRHSElementVect(*projectionSpace.GetFE(i),
+                                     *projectionSpace.GetElementTransformation(i),
+                                     rhs);
+      DenseMatrix rhs_mat(rhs.GetData(), nDofs, projectionSpace.GetVDim());
+      DenseMatrix loc_proj_mat(loc_proj.GetData(), nDofs,
+                               projectionSpace.GetVDim());
+      Mult(invMe, rhs_mat, loc_proj_mat);
+      projectedSol.SetSubVector(dofs, loc_proj);
+   }
+
+   // Compute errors
+   error_estimates.SetSize(nelem);
+   GridFunctionCoefficient proj_gf(&projectedSol);
+   solution.ComputeElementL2Errors(proj_gf, error_estimates);
+   total_error = error_estimates.Norml2();
+}
+
 void LSZienkiewiczZhuEstimator::ComputeEstimates()
 {
-   total_error = LSZZErrorEstimator(integ,
-                                    solution,
-                                    error_estimates,
-                                    subdomain_reconstruction,
-                                    with_coeff,
-                                    tichonov_coeff);
+   total_error =
+      LSZZErrorEstimator(integ, solution, error_estimates,
+                         subdomain_reconstruction, with_coeff, tichonov_coeff);
 
-   current_sequence = solution.FESpace()->GetMesh()->GetSequence();
+   current_sequence = solution.FESpace()->GetSequence();
 }
 
 #ifdef MFEM_USE_MPI
@@ -53,42 +150,44 @@ void L2ZienkiewiczZhuEstimator::ComputeEstimates()
    const double solver_tol = 1e-12;
    const int solver_max_it = 200;
    total_error = L2ZZErrorEstimator(integ, solution, *smooth_flux_space,
-                                    *flux_space, error_estimates,
-                                    local_norm_p, solver_tol, solver_max_it);
+                                    *flux_space, error_estimates, local_norm_p,
+                                    solver_tol, solver_max_it);
 
-   current_sequence = solution.FESpace()->GetMesh()->GetSequence();
+   current_sequence = solution.FESpace()->GetSequence();
 }
 
-#endif // MFEM_USE_MPI
+#endif  // MFEM_USE_MPI
 
-KellyErrorEstimator::KellyErrorEstimator(BilinearFormIntegrator& di_,
-                                         GridFunction& sol_,
-                                         FiniteElementSpace& flux_fespace_,
+KellyErrorEstimator::KellyErrorEstimator(BilinearFormIntegrator &di_,
+                                         GridFunction &sol_,
+                                         FiniteElementSpace &flux_fespace_,
                                          const Array<int> &attributes_)
-   : attributes(attributes_)
-   , flux_integrator(&di_)
-   , solution(&sol_)
-   , flux_space(&flux_fespace_)
-   , own_flux_fespace(false)
+   : attributes(attributes_),
+     flux_integrator(&di_),
+     solution(&sol_),
+     flux_space(&flux_fespace_),
+     own_flux_fespace(false)
 #ifdef MFEM_USE_MPI
-   , isParallel(dynamic_cast<ParFiniteElementSpace*>(sol_.FESpace()))
-#endif // MFEM_USE_MPI
+   ,
+     isParallel(dynamic_cast<ParFiniteElementSpace *>(sol_.FESpace()))
+#endif  // MFEM_USE_MPI
 {
    ResetCoefficientFunctions();
 }
 
-KellyErrorEstimator::KellyErrorEstimator(BilinearFormIntegrator& di_,
-                                         GridFunction& sol_,
-                                         FiniteElementSpace* flux_fespace_,
+KellyErrorEstimator::KellyErrorEstimator(BilinearFormIntegrator &di_,
+                                         GridFunction &sol_,
+                                         FiniteElementSpace *flux_fespace_,
                                          const Array<int> &attributes_)
-   : attributes(attributes_)
-   , flux_integrator(&di_)
-   , solution(&sol_)
-   , flux_space(flux_fespace_)
-   , own_flux_fespace(true)
+   : attributes(attributes_),
+     flux_integrator(&di_),
+     solution(&sol_),
+     flux_space(flux_fespace_),
+     own_flux_fespace(true)
 #ifdef MFEM_USE_MPI
-   , isParallel(dynamic_cast<ParFiniteElementSpace*>(sol_.FESpace()))
-#endif // MFEM_USE_MPI
+   ,
+     isParallel(dynamic_cast<ParFiniteElementSpace *>(sol_.FESpace()))
+#endif  // MFEM_USE_MPI
 {
    ResetCoefficientFunctions();
 }
@@ -103,12 +202,9 @@ KellyErrorEstimator::~KellyErrorEstimator()
 
 void KellyErrorEstimator::ResetCoefficientFunctions()
 {
-   compute_element_coefficient = [](Mesh* mesh, const int e)
-   {
-      return 1.0;
-   };
+   compute_element_coefficient = [](Mesh *mesh, const int e) { return 1.0; };
 
-   compute_face_coefficient = [](Mesh* mesh, const int f,
+   compute_face_coefficient = [](Mesh *mesh, const int f,
                                  const bool shared_face)
    {
       auto FT = [&]()
@@ -116,9 +212,9 @@ void KellyErrorEstimator::ResetCoefficientFunctions()
 #ifdef MFEM_USE_MPI
          if (shared_face)
          {
-            return dynamic_cast<ParMesh*>(mesh)->GetSharedFaceTransformations(f);
+            return dynamic_cast<ParMesh *>(mesh)->GetSharedFaceTransformations(f);
          }
-#endif // MFEM_USE_MPI
+#endif  // MFEM_USE_MPI
          return mesh->GetFaceElementTransformations(f);
       }();
       const auto order = FT->GetFE()->GetOrder();
@@ -148,7 +244,7 @@ void KellyErrorEstimator::ResetCoefficientFunctions()
             diameter = std::max<double>(diameter, p2.DistanceTo(p1));
          }
       }
-      return diameter/(2.0*order);
+      return diameter / (2.0 * order);
    };
 }
 
@@ -179,9 +275,10 @@ void KellyErrorEstimator::ComputeEstimates()
    // 1. Compute fluxes in discontinuous space
    GridFunction *flux =
 #ifdef MFEM_USE_MPI
-      isParallel ? new ParGridFunction(dynamic_cast<ParFiniteElementSpace*>
-                                       (flux_space)) :
-#endif // MFEM_USE_MPI
+      isParallel ? new ParGridFunction(
+         dynamic_cast<ParFiniteElementSpace *>(flux_space))
+      :
+#endif  // MFEM_USE_MPI
       new GridFunction(flux_space);
 
    *flux = 0.0;
@@ -205,7 +302,7 @@ void KellyErrorEstimator::ComputeEstimates()
       xfes->GetElementVDofs(e, xdofs);
       solution->GetSubVector(xdofs, el_x);
 
-      ElementTransformation* Transf = xfes->GetElementTransformation(e);
+      ElementTransformation *Transf = xfes->GetElementTransformation(e);
       flux_integrator->ComputeElementFlux(*xfes->GetFE(e), *Transf, el_x,
                                           *flux_space->GetFE(e), el_f, true);
 
@@ -218,26 +315,34 @@ void KellyErrorEstimator::ComputeEstimates()
    {
       auto FT = mesh->GetFaceElementTransformations(f);
 
-      auto &int_rule = IntRules.Get(FT->FaceGeom, 2 * xfes->GetFaceOrder(f));
-      const auto nip = int_rule.GetNPoints();
+      // auto &int_rule = IntRules.Get(FT->FaceGeom, 2 * xfes->GetFaceOrder(f));
+      // const auto nip = int_rule.GetNPoints();
 
       if (mesh->FaceIsInterior(f))
       {
          int Inf1, Inf2, NCFace;
          mesh->GetFaceInfos(f, &Inf1, &Inf2, &NCFace);
+         int el1, el2;
+         mesh->GetFaceElements(f, &el1, &el2);
+         auto &int_rule =
+            IntRules.Get(FT->FaceGeom,
+                         xfes->GetElementOrder(el1) + xfes->GetElementOrder(el2));
+         const auto nip = int_rule.GetNPoints();
+         // auto &int_rule = IntRules.Get(FT->FaceGeom,
+         // xfes->GetElementOrder(Inf1))
 
          // Convention
          // * Conforming face: Face side with smaller element id handles
          // the integration
          // * Non-conforming face: The slave handles the integration.
          // See FaceInfo documentation for details.
-         bool isNCSlave    = FT->Elem2No >= 0 && NCFace >= 0;
+         bool isNCSlave = FT->Elem2No >= 0 && NCFace >= 0;
          bool isConforming = FT->Elem2No >= 0 && NCFace == -1;
          if ((FT->Elem1No < FT->Elem2No && isConforming) || isNCSlave)
          {
             if (attributes.Size() &&
-                (attributes.FindSorted(FT->Elem1->Attribute) == -1
-                 || attributes.FindSorted(FT->Elem2->Attribute) == -1))
+                (attributes.FindSorted(FT->Elem1->Attribute) == -1 ||
+                 attributes.FindSorted(FT->Elem2->Attribute) == -1))
             {
                continue;
             }
@@ -314,7 +419,7 @@ void KellyErrorEstimator::ComputeEstimates()
                jumps(i) *= jumps(i);
             }
             auto h_k_face = compute_face_coefficient(mesh, f, false);
-            double jump_integral = h_k_face*jumps.Sum();
+            double jump_integral = h_k_face * jumps.Sum();
 
             // A local face is shared between two local elements, so we
             // can get away with integrating the jump only once and add
@@ -326,11 +431,11 @@ void KellyErrorEstimator::ComputeEstimates()
       }
    }
 
-   current_sequence = solution->FESpace()->GetMesh()->GetSequence();
+   current_sequence = solution->FESpace()->GetSequence();
 
 #ifdef MFEM_USE_MPI
    if (!isParallel)
-#endif // MFEM_USE_MPI
+#endif  // MFEM_USE_MPI
    {
       // Finalize element errors
       for (int e = 0; e < xfes->GetNE(); e++)
@@ -350,10 +455,10 @@ void KellyErrorEstimator::ComputeEstimates()
    // 3. Add error contribution from shared interior faces
    // Synchronize face data.
 
-   ParGridFunction *pflux = dynamic_cast<ParGridFunction*>(flux);
+   ParGridFunction *pflux = dynamic_cast<ParGridFunction *>(flux);
    MFEM_VERIFY(pflux, "flux is not a ParGridFunction pointer");
 
-   ParMesh *pmesh = dynamic_cast<ParMesh*>(mesh);
+   ParMesh *pmesh = dynamic_cast<ParMesh *>(mesh);
    MFEM_VERIFY(pmesh, "mesh is not a ParMesh pointer");
 
    pflux->ExchangeFaceNbrData();
@@ -362,8 +467,8 @@ void KellyErrorEstimator::ComputeEstimates()
    {
       auto FT = pmesh->GetSharedFaceTransformations(sf, true);
       if (attributes.Size() &&
-          (attributes.FindSorted(FT->Elem1->Attribute) == -1
-           || attributes.FindSorted(FT->Elem2->Attribute) == -1))
+          (attributes.FindSorted(FT->Elem1->Attribute) == -1 ||
+           attributes.FindSorted(FT->Elem2->Attribute) == -1))
       {
          continue;
       }
@@ -442,7 +547,7 @@ void KellyErrorEstimator::ComputeEstimates()
          jumps(i) *= jumps(i);
       }
       auto h_k_face = compute_face_coefficient(mesh, sf, true);
-      double jump_integral = h_k_face*jumps.Sum();
+      double jump_integral = h_k_face * jumps.Sum();
 
       error_estimates(FT->Elem1No) += jump_integral;
       // We skip "error_estimates(FT->Elem2No) += jump_integral"
@@ -460,14 +565,14 @@ void KellyErrorEstimator::ComputeEstimates()
    }
 
    // Finish by computing the global error.
-   auto pfes = dynamic_cast<ParFiniteElementSpace*>(xfes);
+   auto pfes = dynamic_cast<ParFiniteElementSpace *>(xfes);
    MFEM_VERIFY(pfes, "xfes is not a ParFiniteElementSpace pointer");
 
-   double process_local_error = pow(error_estimates.Norml2(),2.0);
-   MPI_Allreduce(&process_local_error, &total_error, 1, MPI_DOUBLE,
-                 MPI_SUM, pfes->GetComm());
+   double process_local_error = pow(error_estimates.Norml2(), 2.0);
+   MPI_Allreduce(&process_local_error, &total_error, 1, MPI_DOUBLE, MPI_SUM,
+                 pfes->GetComm());
    total_error = sqrt(total_error);
-#endif // MFEM_USE_MPI
+#endif  // MFEM_USE_MPI
 }
 
 void LpErrorEstimator::ComputeEstimates()
@@ -486,16 +591,35 @@ void LpErrorEstimator::ComputeEstimates()
    }
 #ifdef MFEM_USE_MPI
    total_error = error_estimates.Sum();
-   auto pfes = dynamic_cast<ParFiniteElementSpace*>(sol->FESpace());
+   auto pfes = dynamic_cast<ParFiniteElementSpace *>(sol->FESpace());
    if (pfes)
    {
       auto process_local_error = total_error;
-      MPI_Allreduce(&process_local_error, &total_error, 1, MPI_DOUBLE,
-                    MPI_SUM, pfes->GetComm());
+      MPI_Allreduce(&process_local_error, &total_error, 1, MPI_DOUBLE, MPI_SUM,
+                    pfes->GetComm());
    }
-#endif // MFEM_USE_MPI
-   total_error = pow(total_error, 1.0/local_norm_p);
-   current_sequence = sol->FESpace()->GetMesh()->GetSequence();
+#endif  // MFEM_USE_MPI
+   total_error = pow(total_error, 1.0 / local_norm_p);
+   current_sequence = sol->FESpace()->GetSequence();
 }
 
-} // namespace mfem
+SolJumpEstimator::SolJumpEstimator(GridFunction &sol_)
+   : solution(&sol_) {}
+
+void SolJumpEstimator::ComputeEstimates()
+{
+   const int nelem = solution->FESpace()->GetNE();
+   error_estimates.SetSize(nelem);
+
+   FiniteElementSpace *fespace = solution->FESpace();
+   Mesh *mesh = fespace->GetMesh();
+   FiniteElementSpace *h1fes = new FiniteElementSpace(mesh, new H1_FECollection(fespace->GetMaxElementOrder(), mesh->Dimension()));
+   GridFunctionCoefficient sol_cf(solution);
+   GridFunction h1sol(h1fes);
+   h1sol.ProjectDiscCoefficient(sol_cf, mfem::GridFunction::AvgType::ARITHMETIC);
+   h1sol.ComputeElementL2Errors(sol_cf, error_estimates);
+   total_error = error_estimates.Norml2();
+   delete h1fes;
+}
+
+}  // namespace mfem

fem/estimators.hpp

@@ -70,6 +70,113 @@ public:
    virtual const Array<int> &GetAnisotropicFlags() = 0;
 };
 
+class PRefDiffEstimator : public ErrorEstimator
+{
+public:
+
+private:
+   int current_sequence = -1;
+
+   Vector error_estimates;
+
+   double total_error = 0.0;
+
+   int p_comp = -1;
+
+   GridFunction* solution;               ///< Not owned.
+
+   /// Check if the mesh of the solution was modified.
+   bool MeshIsModified()
+   {
+      long mesh_sequence = solution->FESpace()->GetMesh()->GetSequence();
+      MFEM_ASSERT(mesh_sequence >= current_sequence,
+                  "improper mesh update sequence");
+      return (mesh_sequence > current_sequence);
+   }
+
+   void ComputeEstimates();
+
+public:
+   /** @brief Construct a new PRefDiffEstimator object for a scalar field.
+       @param sol_        The solution field whose error is to be estimated.
+       @param p_comp_     Comparison order.. if -ve, subtract this from current order
+   */
+   PRefDiffEstimator(GridFunction& sol_, int p_comp_);
+
+   ~PRefDiffEstimator() { };
+
+   /// Get a Vector with all element errors.
+   const Vector& GetLocalErrors() override
+   {
+      if (MeshIsModified())
+      {
+         ComputeEstimates();
+      }
+      return error_estimates;
+   }
+
+   /// Reset the error estimator.
+   void Reset() override { current_sequence = -1; };
+
+   virtual double GetTotalError() const override { return total_error; }
+};
+
+/**
+ * @brief Highest
+ *
+ */
+class ProjectionErrorEstimator : public ErrorEstimator
+{
+protected:
+   long current_sequence; // current FE sequence
+   Vector error_estimates; // element-wise error
+   double total_error; // total error
+   int offset; // degree offset
+
+   GridFunction &solution;
+
+   /// Check if the finite element of the solution was modified.
+   bool FESpaceIsModified()
+   {
+      long fe_sequence = solution.FESpace()->GetSequence();
+      MFEM_ASSERT(fe_sequence >= current_sequence, "");
+      return (fe_sequence > current_sequence);
+   }
+
+   /// Compute the element error estimates.
+   void ComputeEstimates();
+
+public:
+   /** @brief Construct a new ProjectionErrorEstimator object. ||π_{k-s}u_h - u_h||
+       @param sol      The solution field whose error is to be estimated.
+       @param degree_offset The degree offset (s) for the projected space. Default=1
+   */
+   ProjectionErrorEstimator(GridFunction &sol, const int degree_offset=1)
+      : current_sequence(-1),
+        total_error(-1.0),
+        offset(degree_offset),
+        solution(sol)
+   { }
+
+   /// Return the total error from the last error estimate.
+   virtual double GetTotalError() const override { return total_error; }
+
+   /// Get a Vector with all element errors.
+   virtual const Vector &GetLocalErrors() override
+   {
+      if (FESpaceIsModified())
+      {
+         ComputeEstimates();
+      }
+      return error_estimates;
+   }
+
+   /// Reset the error estimator.
+   virtual void Reset() override { current_sequence = -1; }
+
+   virtual ~ProjectionErrorEstimator() { }
+};
+
 
 /** @brief The ZienkiewiczZhuEstimator class implements the Zienkiewicz-Zhu
     error estimation procedure.
@@ -104,11 +211,11 @@ protected:
    bool own_flux_fes; ///< Ownership flag for flux_space.
 
    /// Check if the mesh of the solution was modified.
-   bool MeshIsModified()
+   bool FESpaceIsModified()
    {
-      long mesh_sequence = solution.FESpace()->GetMesh()->GetSequence();
-      MFEM_ASSERT(mesh_sequence >= current_sequence, "");
-      return (mesh_sequence > current_sequence);
+      long fe_sequence = solution.FESpace()->GetSequence();
+      MFEM_ASSERT(fe_sequence >= current_sequence, "");
+      return (fe_sequence > current_sequence);
    }
 
    /// Compute the element error estimates.
@@ -177,7 +284,7 @@ public:
    /// Get a Vector with all element errors.
    virtual const Vector &GetLocalErrors() override
    {
-      if (MeshIsModified()) { ComputeEstimates(); }
+      if (FESpaceIsModified()) { ComputeEstimates(); }
       return error_estimates;
    }
 
@@ -186,7 +293,7 @@ public:
        enabled. */
    virtual const Array<int> &GetAnisotropicFlags() override
    {
-      if (MeshIsModified()) { ComputeEstimates(); }
+      if (FESpaceIsModified()) { ComputeEstimates(); }
       return aniso_flags;
    }
 
@@ -252,11 +359,11 @@ protected:
    bool with_coeff;
 
    /// Check if the mesh of the solution was modified.
-   bool MeshIsModified()
+   bool FESpaceIsModified()
    {
-      long mesh_sequence = solution.FESpace()->GetMesh()->GetSequence();
-      MFEM_ASSERT(mesh_sequence >= current_sequence, "");
-      return (mesh_sequence > current_sequence);
+      long fe_sequence = solution.FESpace()->GetSequence();
+      MFEM_ASSERT(fe_sequence >= current_sequence, "");
+      return (fe_sequence > current_sequence);
    }
 
    /// Compute the element error estimates.
@@ -303,7 +410,7 @@ public:
    /// Get a Vector with all element errors.
    virtual const Vector &GetLocalErrors() override
    {
-      if (MeshIsModified()) { ComputeEstimates(); }
+      if (FESpaceIsModified()) { ComputeEstimates(); }
       return error_estimates;
    }
 
@@ -345,11 +452,11 @@ protected:
    bool own_flux_fes; ///< Ownership flag for flux_space and smooth_flux_space.
 
    /// Check if the mesh of the solution was modified.
-   bool MeshIsModified()
+   bool FESpaceIsModified()
    {
-      long mesh_sequence = solution.FESpace()->GetMesh()->GetSequence();
-      MFEM_ASSERT(mesh_sequence >= current_sequence, "");
-      return (mesh_sequence > current_sequence);
+      long fe_sequence = solution.FESpace()->GetSequence();
+      MFEM_ASSERT(fe_sequence >= current_sequence, "");
+      return (fe_sequence > current_sequence);
    }
 
    /// Compute the element error estimates.
@@ -416,7 +523,7 @@ public:
    /// Get a Vector with all element errors.
    virtual const Vector &GetLocalErrors() override
    {
-      if (MeshIsModified()) { ComputeEstimates(); }
+      if (FESpaceIsModified()) { ComputeEstimates(); }
       return error_estimates;
    }
 
@@ -459,11 +566,11 @@ protected:
    GridFunction * sol;
 
    /// Check if the mesh of the solution was modified.
-   bool MeshIsModified()
+   bool FESpaceIsModified()
    {
-      long mesh_sequence = sol->FESpace()->GetMesh()->GetSequence();
-      MFEM_ASSERT(mesh_sequence >= current_sequence, "");
-      return (mesh_sequence > current_sequence);
+      long fe_sequence = sol->FESpace()->GetSequence();
+      MFEM_ASSERT(fe_sequence >= current_sequence, "");
+      return (fe_sequence > current_sequence);
    }
 
    /// Compute the element error estimates.
@@ -510,7 +617,7 @@ public:
    /// Get a Vector with all element errors.
    virtual const Vector &GetLocalErrors() override
    {
-      if (MeshIsModified()) { ComputeEstimates(); }
+      if (FESpaceIsModified()) { ComputeEstimates(); }
       return error_estimates;
    }
 
@@ -603,12 +710,12 @@ private:
 #endif
 
    /// Check if the mesh of the solution was modified.
-   bool MeshIsModified()
+   bool FESpaceIsModified()
    {
-      long mesh_sequence = solution->FESpace()->GetMesh()->GetSequence();
-      MFEM_ASSERT(mesh_sequence >= current_sequence,
+      long fe_sequence = solution->FESpace()->GetSequence();
+      MFEM_ASSERT(fe_sequence >= current_sequence,
                   "improper mesh update sequence");
-      return (mesh_sequence > current_sequence);
+      return (fe_sequence > current_sequence);
    }
 
    /** @brief Compute the element error estimates.
@@ -651,7 +758,7 @@ public:
    /// Get a Vector with all element errors.
    const Vector& GetLocalErrors() override
    {
-      if (MeshIsModified())
+      if (FESpaceIsModified())
       {
          ComputeEstimates();
       }
@@ -690,6 +797,56 @@ public:
    void ResetCoefficientFunctions();
 };
 
+class SolJumpEstimator : public ErrorEstimator
+{
+public:
+
+private:
+   int current_sequence = -1;
+
+   Vector error_estimates;
+
+   double total_error = 0.0;
+
+   GridFunction* solution;               ///< Not owned.
+
+   /// Check if the mesh of the solution was modified.
+   bool MeshIsModified()
+   {
+      long mesh_sequence = solution->FESpace()->GetMesh()->GetSequence();
+      MFEM_ASSERT(mesh_sequence >= current_sequence,
+                  "improper mesh update sequence");
+      return (mesh_sequence > current_sequence);
+   }
+
+   void ComputeEstimates();
+
+public:
+   /** @brief Construct a new SolJumpEstimator object for a scalar field.
+       @param sol_        The solution field whose error is to be estimated.
+       @param p_comp_     Comparison order.. if -ve, subtract this from current order
+   */
+   SolJumpEstimator(GridFunction& sol_);
+
+   ~SolJumpEstimator() { };
+
+   /// Get a Vector with all element errors.
+   const Vector& GetLocalErrors() override
+   {
+      if (MeshIsModified())
+      {
+         ComputeEstimates();
+      }
+      return error_estimates;
+   }
+
+   /// Reset the error estimator.
+   void Reset() override { current_sequence = -1; };
+
+   virtual double GetTotalError() const override { return total_error; }
+};
+
+
 } // namespace mfem
 
 #endif // MFEM_ERROR_ESTIMATORS

fem/gridfunc.cpp

@@ -4583,4 +4583,47 @@ GridFunction *Extrude1DGridFunction(Mesh *mesh, Mesh *mesh2d,
    return sol2d;
 }
 
+
+GridFunction* ProlongToMaxOrder(const GridFunction *x, const int fieldtype)
+{
+   const FiniteElementSpace *fespace = x->FESpace();
+   Mesh *mesh = fespace->GetMesh();
+   const FiniteElementCollection *fec = fespace->FEColl();
+
+   // find the max order in the space
+   const int max_order = fespace->GetMaxElementOrder();
+
+   // create a visualization space of max order for all elements
+   FiniteElementCollection *fecInt = NULL;
+   if (fieldtype == 0)
+   {
+      fecInt = new H1_FECollection(max_order, mesh->Dimension());
+   }
+   else if (fieldtype == 1)
+   {
+      fecInt = new L2_FECollection(max_order, mesh->Dimension());
+   }
+   FiniteElementSpace *spaceInt = new FiniteElementSpace(mesh, fecInt,
+                                                         fespace->GetVDim());
+
+   IsoparametricTransformation T;
+   DenseMatrix I;
+
+   GridFunction *xInt = new GridFunction(spaceInt);
+
+   if (fespace->GetVDim() == 1)
+   {
+      GridFunctionCoefficient cf(x);
+      xInt->ProjectCoefficient(cf);
+   }
+   else
+   {
+      VectorGridFunctionCoefficient cf(x);
+      xInt->ProjectCoefficient(cf);
+   }
+
+   xInt->MakeOwner(fecInt);
+   return xInt;
+}
+
 }

fem/gridfunc.hpp

@@ -877,6 +877,15 @@ public:
 GridFunction *Extrude1DGridFunction(Mesh *mesh, Mesh *mesh2d,
                                     GridFunction *sol, const int ny);
 
+/**
+ * @brief Prolongate varying-order GridFunction to equal-order GridFunction by interpolation. Usually used for visualization purpose.
+ *
+ * @param x Varying order GridFunction
+ * @return GridFunction Max-order GridFunction
+ */
+// Experimental - required for visualizing functions on p-refined spaces.
+GridFunction* ProlongToMaxOrder(const GridFunction *x, const int fieldtype=1);
+
 } // namespace mfem
 
 #endif

fem/hyperbolic_conservation_laws.cpp

@@ -0,0 +1,39 @@
+#include "hyperbolic_conservation_laws.hpp"
+
+namespace mfem
+{
+void FluxReconErrorEstimator::ComputeEstimates()
+{
+   // FiniteElementSpace* fespace = solution.FESpace();
+   // Mesh* mesh = fespace->GetMesh();
+   // const int sdim = mesh->Dimension();
+   // const int vdim = fespace->GetVDim();
+   // const int numElem = mesh->GetNE();
+   // const int numFace = mesh->GetNumFaces();
+   // const int base_order = fespace->GetElementOrder(0);
+   // FiniteElementSpace RT(mesh, RTfec);
+   // for (int i=0; i<numElem; i++)
+   // {
+   //    RT.SetElementOrder(i, fespace->GetElementOrder(i));
+   // }
+   // RT.Update(false);
+
+   // // Create gridfunction for each row of the flux
+   // Array<GridFunction*> fluxes(vdim);
+   // for (auto &flux:fluxes)
+   // {
+   //    flux = new GridFunction(&RT);
+   // }
+
+   // Array<int> dofs;
+   // DenseMatrix flux_val;
+   // for (int i=0; i< numFace; i++)
+   // {
+   //    FiniteElement* fe = RT->GetFE(i);
+      
+
+   // }
+   // delete fluxes;
+   mfem_error("WIP");
+}
+}

fem/proxGalerkinHCL.hpp

@@ -0,0 +1,531 @@
+#pragma once
+
+#include <fstream>
+#include <iostream>
+#include <sstream>
+#include <cmath>
+#include <functional>
+
+#include "mfem.hpp"
+#include "fem/hyperbolic_conservation_laws.hpp"
+
+namespace mfem
+{
+class MFEMNew
+{
+public:
+   static NonlinearForm* newNonlinearForm(FiniteElementSpace *fes)
+   {
+#ifdef MFEM_USE_MPI
+      auto *pfes = dynamic_cast<ParFiniteElementSpace*>(fes);
+      if (pfes)
+      {
+         return new ParNonlinearForm(pfes);
+      }
+      else
+      {
+         return new NonlinearForm(fes);
+      }
+#else
+      return new NonlinearForm(fes);
+#endif
+   }
+   static BilinearForm* newBilinearForm(FiniteElementSpace* fes)
+   {
+#ifdef MFEM_USE_MPI
+      auto *pfes = dynamic_cast<ParFiniteElementSpace*>(fes);
+      if (pfes)
+      {
+         return new ParBilinearForm(pfes);
+      }
+      else
+      {
+         return new BilinearForm(fes);
+      }
+#else
+      return new BilinearForm(fes);
+#endif
+   }
+   static LinearForm* newLinearForm(FiniteElementSpace* fes)
+   {
+#ifdef MFEM_USE_MPI
+      auto *pfes = dynamic_cast<ParFiniteElementSpace*>(fes);
+      if (pfes)
+      {
+         return new ParLinearForm(pfes);
+      }
+      else
+      {
+         return new LinearForm(fes);
+      }
+#else
+      return new LinearForm(fes);
+#endif
+   }
+   static LinearForm* newLinearForm(FiniteElementSpace* fes, double *data)
+   {
+#ifdef MFEM_USE_MPI
+      auto *pfes = dynamic_cast<ParFiniteElementSpace*>(fes);
+      if (pfes)
+      {
+         return new ParLinearForm(pfes, data);
+      }
+      else
+      {
+         return new LinearForm(fes, data);
+      }
+#else
+      return new LinearForm(fes, data);
+#endif
+   }
+   static GridFunction* newGridFunction(FiniteElementSpace* fes)
+   {
+#ifdef MFEM_USE_MPI
+      auto *pfes = dynamic_cast<ParFiniteElementSpace*>(fes);
+      if (pfes)
+      {
+         return new ParGridFunction(pfes);
+      }
+      else
+      {
+         return new GridFunction(fes);
+      }
+#else
+      return new GridFunction(fes);
+#endif
+   }
+   static GridFunction* newGridFunction(FiniteElementSpace* fes, double *data)
+   {
+#ifdef MFEM_USE_MPI
+      auto *pfes = dynamic_cast<ParFiniteElementSpace*>(fes);
+      if (pfes)
+      {
+         return new ParGridFunction(pfes, data);
+      }
+      else
+      {
+         return new GridFunction(fes, data);
+      }
+#else
+      return new GridFunction(fes, data);
+#endif
+   }
+   static GridFunction* newGridFunction(FiniteElementSpace* fes, Vector &base,
+                                        int offset=0)
+   {
+#ifdef MFEM_USE_MPI
+      auto *pfes = dynamic_cast<ParFiniteElementSpace*>(fes);
+      if (pfes)
+      {
+         return new ParGridFunction(pfes, base, offset);
+      }
+      else
+      {
+         return new GridFunction(fes, base, offset);
+      }
+#else
+      return new GridFunction(fes, base, offset);
+#endif
+   }
+};
+class Scalar2ScalarMappedGF : public Coefficient
+{
+public:
+   Scalar2ScalarMappedGF(GridFunction *gf,
+                         std::function<double(double, ElementTransformation &, const IntegrationPoint &)>
+                         fun,
+                         int comp=1)
+      :gfc(gf, comp), fun(fun) {}
+   virtual double Eval(ElementTransformation &T,
+                       const IntegrationPoint &ip)
+   {
+      return fun(gfc.Eval(T, ip), T, ip);
+   }
+protected:
+private:
+   GridFunctionCoefficient gfc;
+   std::function<double(double, ElementTransformation &, const IntegrationPoint &)>
+   fun;
+};
+class Vector2ScalarMappedGF : public Coefficient
+{
+public:
+   Vector2ScalarMappedGF(GridFunction *gf,
+                         std::function<double(Vector&, ElementTransformation &, const IntegrationPoint &)>
+                         fun)
+      :gfc(gf), fun(fun), gf_vdim(gf->VectorDim())
+   {
+#ifndef MFEM_THREAD_SAFE
+      val.SetSize(gf_vdim);
+#endif
+   }
+   virtual double Eval(ElementTransformation &T,
+                       const IntegrationPoint &ip)
+   {
+#ifdef MFEM_THREAD_SAFE
+      Vector val(gf_vdim);
+#endif
+      gfc.Eval(val, T, ip);
+      return fun(val, T, ip);
+   }
+protected:
+private:
+#ifndef MFEM_THREAD_SAFE
+   Vector val;
+#endif
+   VectorGridFunctionCoefficient gfc;
+   const int gf_vdim;
+   std::function<double(Vector&, ElementTransformation &, const IntegrationPoint &)>
+   fun;
+};
+
+
+
+class Scalar2VectorMappedGF : public VectorCoefficient
+{
+public:
+   Scalar2VectorMappedGF(GridFunction *gf,
+                         std::function<void(Vector&, double, ElementTransformation &, const IntegrationPoint &)>
+                         fun, int vd, int comp=1)
+      :VectorCoefficient(vd), gfc(gf, comp), fun(fun) {}
+   virtual void Eval(Vector &V, ElementTransformation &T,
+                     const IntegrationPoint &ip)
+   {
+      V.SetSize(vdim);
+      fun(V, gfc.Eval(T, ip), T, ip);
+      return;
+   }
+protected:
+private:
+   GridFunctionCoefficient gfc;
+   std::function<void(Vector&, double, ElementTransformation &, const IntegrationPoint &)>
+   fun;
+};
+class Vector2VectorMappedGF : public VectorCoefficient
+{
+public:
+   Vector2VectorMappedGF(GridFunction *gf,
+                         std::function<void(Vector&, const Vector&, ElementTransformation &, const IntegrationPoint &)>
+                         fun, int vd)
+      :VectorCoefficient(vd), gfc(gf), fun(fun), gf_vdim(gf->VectorDim())
+   {
+#ifndef MFEM_THREAD_SAFE
+      val.SetSize(gf_vdim);
+#endif
+   }
+   virtual void Eval(Vector &V, ElementTransformation &T,
+                     const IntegrationPoint &ip)
+   {
+      V.SetSize(vdim);
+#ifdef MFEM_THREAD_SAFE
+      Vector val(gf_vdim);
+#endif
+      gfc.Eval(val, T, ip);
+      fun(V, val, T, ip);
+      return;
+   }
+protected:
+private:
+#ifndef MFEM_THREAD_SAFE
+   Vector val;
+#endif
+   VectorGridFunctionCoefficient gfc;
+   const int gf_vdim;
+   std::function<void(Vector&, const Vector&, ElementTransformation &, const IntegrationPoint &)>
+   fun;
+};
+
+class Scalar2MatrixMappedGF : public MatrixCoefficient
+{
+public:
+   Scalar2MatrixMappedGF(GridFunction *gf,
+                         std::function<void(DenseMatrix&, double)> fun, int h, int w, int comp=1)
+      :MatrixCoefficient(h,w), gfc(gf, comp), fun(fun) {}
+   virtual void Eval(DenseMatrix &M, ElementTransformation &T,
+                     const IntegrationPoint &ip)
+   {
+      M.SetSize(height, width);
+      fun(M, gfc.Eval(T, ip));
+      return;
+   }
+protected:
+private:
+   GridFunctionCoefficient gfc;
+   std::function<void(DenseMatrix&, double)> fun;
+};
+class Vector2MatrixMappedGF : public MatrixCoefficient
+{
+public:
+   Vector2MatrixMappedGF(GridFunction *gf,
+                         std::function<void(DenseMatrix&, const Vector&, ElementTransformation &, const IntegrationPoint &)>
+                         fun, int h, int w)
+      :MatrixCoefficient(h, w), gfc(gf), fun(fun), gf_vdim(gf->VectorDim())
+   {
+#ifndef MFEM_THREAD_SAFE
+      val.SetSize(gf_vdim);
+#endif
+   }
+   virtual void Eval(DenseMatrix &M, ElementTransformation &T,
+                     const IntegrationPoint &ip)
+   {
+      M.SetSize(height, width);
+#ifdef MFEM_THREAD_SAFE
+      Vector val(gf_vdim);
+#endif
+      gfc.Eval(val, T, ip);
+      fun(M, val, T, ip);
+      return;
+   }
+protected:
+private:
+#ifndef MFEM_THREAD_SAFE
+   Vector val;
+#endif
+   VectorGridFunctionCoefficient gfc;
+   const int gf_vdim;
+   std::function<void(DenseMatrix&, const Vector&, ElementTransformation &, const IntegrationPoint &)>
+   fun;
+};
+
+class ProxGalerkinAlphaMaker
+{
+public:
+   ProxGalerkinAlphaMaker() {}
+   virtual double GetAlpha(const int k) {mfem_error("Not implemented"); return 0.0;}
+};
+
+class ProxGalerkinPolynomialAlphaMaker : public ProxGalerkinAlphaMaker
+{
+public:
+   ProxGalerkinPolynomialAlphaMaker(const double alpha0,
+                                    const double degree):alpha0(alpha0), degree(degree),
+      ProxGalerkinAlphaMaker() {}
+   virtual double GetAlpha(const int k) {return alpha0*std::pow(k, degree);}
+protected:
+   const double alpha0;
+   const double degree;
+};
+
+class ProxGalerkinExponentialAlphaMaker : public ProxGalerkinAlphaMaker
+{
+public:
+   ProxGalerkinExponentialAlphaMaker(const double alpha0,
+                                     const double base):alpha0(alpha0), base(base),
+      ProxGalerkinAlphaMaker() {}
+   virtual double GetAlpha(const int k) {return alpha0*std::pow(base, k);}
+protected:
+   const double alpha0;
+   const double base;
+};
+
+class ProxGalerkinHCL : public TimeDependentOperator
+{
+private:
+   // Vector finite element space containing conserved variables
+   FiniteElementSpace *fes;
+   FiniteElementSpace *vfes;
+   // Element integration form. Should contain ComputeFlux
+   HyperbolicFormIntegrator &formIntegrator;
+   // Base Nonlinear Form
+   NonlinearForm *nonlinearForm;
+   // element-wise inverse mass matrix
+   // std::vector<DenseMatrix> Me_inv;
+   BilinearForm *M_inv;
+   // global maximum characteristic speed. Updated by form integrators
+   mutable double max_char_speed;
+   // auxiliary variable used in Mult
+   mutable Vector z;
+   ProxGalerkinAlphaMaker &alphamaker;
+   GridFunction *latent_k, *latent, *delta_latent, *delta_dxdt;
+   BilinearForm *M;
+   NonlinearForm *latentMinv;
+   Vector2VectorMappedGF FLatent;
+   Vector2VectorMappedGF NewtonFLatent;
+   const int dim;
+   const int num_equations;
+   const int maxit = 1e03;
+
+   // Compute element-wise inverse mass matrix
+   void ComputeInvMass();
+
+public:
+   /**
+    * @brief Construct a new ProxGalerkinHCL object
+    *
+    * @param vfes_ vector finite element space. Only tested for DG [Pₚ]ⁿ
+    * @param formIntegrator_ (F(u,x), grad v) and (F̂(u±, x, n), [[v]])
+    * @param num_equations_ the number of equations
+    */
+   ProxGalerkinHCL(
+      FiniteElementSpace *fes,
+      FiniteElementSpace *vfes,
+      HyperbolicFormIntegrator &formIntegrator,
+      const int num_equations,
+      ProxGalerkinAlphaMaker &alphamaker);
+   /**
+    * @brief Apply nonlinear form to obtain M⁻¹(DIVF + JUMP HAT(F))
+    *
+    * @param x current solution vector
+    * @param y resulting dual vector to be used in an EXPLICIT solver
+    */
+   virtual void Mult(const Vector &x, Vector &y) const;
+   // get global maximum characteristic speed to be used in CFL condition
+   // where max_char_speed is updated during Mult.
+   inline double getMaxCharSpeed()
+   {
+      return max_char_speed;
+   }
+
+   void Update()
+   {
+      nonlinearForm->Update();
+      height = nonlinearForm->Height();
+      width = height;
+      z.SetSize(height);
+
+      ComputeInvMass();
+   }
+   virtual void ImplicitSolve(const double dt, const Vector &x, Vector &dxdt);
+
+   virtual ~ProxGalerkinHCL()
+   {
+      delete latent_k; delete latent; delete delta_dxdt; delta_latent;
+      delete nonlinearForm; delete latentMinv; delete M_inv; delete M;
+   }
+};
+
+//////////////////////////////////////////////////////////////////
+///        HYPERBOLIC CONSERVATION LAWS IMPLEMENTATION         ///
+//////////////////////////////////////////////////////////////////
+
+// Implementation of class ProxGalerkinHCL
+ProxGalerkinHCL::ProxGalerkinHCL(
+   FiniteElementSpace *fes,
+   FiniteElementSpace *vfes,
+   HyperbolicFormIntegrator &formIntegrator,
+   const int num_equations,
+   ProxGalerkinAlphaMaker &alphamaker)
+   : fes(fes),
+     vfes(vfes),
+     formIntegrator(formIntegrator),
+     M_inv(nullptr),
+     M(nullptr),
+     z(vfes->GetVSize()),
+     dim(vfes->GetFE(0)->GetDim()),
+     num_equations(num_equations),
+     alphamaker(alphamaker),
+     latent_k(MFEMNew::newGridFunction(vfes)),
+     latent(MFEMNew::newGridFunction(vfes)),
+     delta_dxdt(MFEMNew::newGridFunction(vfes)),
+     delta_latent(MFEMNew::newGridFunction(vfes)),
+     FLatent(latent, [](Vector &y, const Vector &x, ElementTransformation &T,
+                        const IntegrationPoint &ip) {y = std::exp(x[0]); }, latent->VectorDim()),
+NewtonFLatent(latent, [](Vector &y, const Vector &x, ElementTransformation &T,
+const IntegrationPoint &ip) {y = -std::exp(x[0])*(1.0-x[0]); },
+latent->VectorDim()),
+TimeDependentOperator(vfes->GetVSize())
+{
+   // Standard local assembly and inversion for energy mass matrices.
+   ComputeInvMass();
+   nonlinearForm = MFEMNew::newNonlinearForm(vfes);
+   formIntegrator.resetMaxCharSpeed();
+
+   nonlinearForm->AddDomainIntegrator(&formIntegrator);
+   nonlinearForm->AddInteriorFaceIntegrator(&formIntegrator);
+   latentMinv = MFEMNew::newNonlinearForm(vfes);
+   latentMinv->AddDomainIntegrator(new InverseIntegrator(new
+                                                         VectorMassIntegrator(FLatent)));
+}
+
+void ProxGalerkinHCL::ComputeInvMass()
+{
+   if (M_inv) { delete M_inv; }
+   if (M) {delete M; }
+   M_inv = MFEMNew::newBilinearForm(vfes);
+   M_inv->AddDomainIntegrator(new InverseIntegrator(new VectorMassIntegrator()));
+   M_inv->Assemble();
+   M = MFEMNew::newBilinearForm(vfes);
+   M->AddDomainIntegrator(new VectorMassIntegrator());
+   M->Assemble();
+}
+
+void ProxGalerkinHCL::Mult(const Vector &x, Vector &y) const
+{
+   // 0. Reset wavespeed computation before operator application.
+   formIntegrator.resetMaxCharSpeed();
+   // 1. Create the vector z with the face terms (F(u), grad v) - <F.n(u), [w]>.
+   nonlinearForm->Mult(x, z);
+   max_char_speed = formIntegrator.getMaxCharSpeed();
+
+   // 2. Multiply element-wise by the inverse mass matrices.
+   int dof = M_inv->Size();
+   Vector zval(z.GetData(), dof);
+   Vector yval(y.GetData(), dof);
+   for (int i=0; i<num_equations; i++)
+   {
+      zval.SetData(z.GetData() + i*dof);
+      yval.SetData(y.GetData() + i*dof);
+      M_inv->Mult(zval, yval);
+   }
+}
+
+void ProxGalerkinHCL::ImplicitSolve(const double dt, const Vector &x,
+                                    Vector &dxdt)
+{
+
+   *latent = x;
+   latent->ApplyMap([](double x) {return std::log(x);});
+   GridFunction *xnew = MFEMNew::newGridFunction(vfes);
+   GridFunction *Mrhs = MFEMNew::newGridFunction(vfes);
+   LinearForm *newtonRHS = MFEMNew::newLinearForm(vfes);
+   newtonRHS->AddDomainIntegrator(new VectorDomainLFIntegrator(NewtonFLatent));
+   VectorGridFunctionCoefficient xnew_cf(xnew);
+   newtonRHS->AddDomainIntegrator(new VectorDomainLFIntegrator(xnew_cf));
+   *delta_dxdt = infinity();
+   *delta_latent = infinity();
+   dxdt.SetSize(x.Size());
+   dxdt = 0.0;
+   *xnew = x;
+   xnew->Add(dt, dxdt);
+   int k=0;
+   while (k++ < maxit & delta_latent->Normlinf() > 1e-07)
+   {
+      double alpha = alphamaker.GetAlpha(k);
+      *latent_k = *latent;
+      out << "Step " << k << ", " << alpha << std::endl;
+      while (delta_dxdt->Normlinf() > 1e-07)
+      {
+         *delta_dxdt = dxdt;
+
+         Mult(*xnew, dxdt);
+         dxdt.Add(alpha, *latent_k);
+         dxdt.Add(-alpha, *latent);
+
+         *delta_dxdt -= dxdt;
+
+         *xnew = x;
+         xnew->Add(dt, dxdt);
+
+         *delta_latent = infinity();
+         while (delta_latent->Normlinf() > 1e-08)
+         {
+            *delta_latent = *latent;
+            newtonRHS->Assemble();
+            M->Mult(*newtonRHS, *Mrhs);
+            latentMinv->Mult(*Mrhs, *latent);
+            *delta_latent -= *latent;
+            out << latent->Normlinf() << std::endl;
+         }
+
+         out << delta_dxdt->Normlinf() << ", " << delta_latent->Normlinf() << std::endl;
+      }
+      *delta_latent = *latent_k;
+      *delta_latent -= *latent;
+   }
+
+   delete xnew;
+   delete Mrhs;
+   delete newtonRHS;
+}
+
+} // namespace mfem

linalg/densemat.cpp

@@ -4211,6 +4211,7 @@ void DenseMatrixSVD::Eval(DenseMatrix &M)
 
    if (info)
    {
+      M.Print();
       mfem::err << "DenseMatrixSVD::Eval() : info = " << info << endl;
       mfem_error();
    }

linalg/vector.cpp

@@ -312,6 +312,13 @@ void Vector::Neg()
    auto y = ReadWrite(use_dev);
    mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i) { y[i] = -y[i]; });
 }
+void Vector::ApplyMap(std::function<double(double)> fun)
+{
+   const bool use_dev = UseDevice();
+   const int N = size;
+   auto y = ReadWrite(use_dev);
+   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i) { y[i] = fun(y[i]); });
+}
 
 void Vector::Reciprocal()
 {

linalg/vector.hpp

@@ -25,6 +25,7 @@
 #include <cmath>
 #include <iostream>
 #include <limits>
+#include <functional>
 #if defined(_MSC_VER) && (_MSC_VER < 1800)
 #include <float.h>
 #define isfinite _finite
@@ -323,6 +324,8 @@ public:
    /// (*this) = -(*this)
    void Neg();
 
+   void ApplyMap(std::function<double(double)> fun);
+
    /// (*this)(i) = 1.0 / (*this)(i)
    void Reciprocal();