Merge branch 'master' into HCL-refactor-Jacobian

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/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/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");
+}
+}