Fix another makefile bug

miniapps/conslaw/.gitignore

@@ -0,0 +1,4 @@
+build
+lib
+*.mesh
+*.gf

miniapps/conslaw/README.md

@@ -0,0 +1,35 @@
+# Discontinuous Galerkin MFEM mini-app using partial assembly
+
+This mini-app demonstrates using partial-assembly to solve hyperbolic
+conservation laws using discontinuous Galerkin methods and explicit
+time integration.
+
+The main object is a `PartialAssembly` object, which provides:
+
+- local interpolation and differentiation operators (including at faces)
+- face access to metric terms at quadrature data
+
+On top of this object, there are several operators that are provided
+(but more are possible). All of these operators allow for a "coefficient"
+to be evaluated at quadrature points, which is referred to as D. These
+operators are:
+
+- `BtDB`, which represents mass or source terms with coefficient `D`
+- `GtDB`, which represents dot product with the gradient of test functions
+- `BtDB_face`, which represents integrating against test functions on faces
+
+Using any of these operators simply requires templating on a class `D` which
+provides an operator to evaluate the coefficient at a quadrature point.
+
+There is a `ConservationLaw` object which is built on these three operators.
+Given a flux function and numerical flux function, it will assemble the 
+corresponding DG residual.
+
+Examples are provided in the `apps` directory for solving the scalar advection
+equation, Burgers' equation and the Euler equations of gas dynamics.
+
+This mini-app is still incomplete. Improvements are needed for:
+
+[ ] Handing of mixed meshes
+[ ] AMR and non-conforming meshes
+[ ] Second-order operators and viscous terms

miniapps/conslaw/apps/.gitignore

@@ -0,0 +1,3 @@
+*
+!*/
+!*.*

miniapps/conslaw/apps/advection.cpp

@@ -0,0 +1,208 @@
+//          MFEM DG FOR CONSERVATION LAWS WITH PARTIAL ASSEMBLY
+
+#include "mfem.hpp"
+#include "dg.hpp"
+#include <fstream>
+#include <iostream>
+#include <random>
+
+#include "advection.hpp"
+#include "evol.hpp"
+
+using namespace std;
+using namespace mfem;
+using namespace dg;
+
+
+Vector bb_min, bb_max;
+double u0_function(const Vector &x);
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   const char *mesh_file = "../../data/periodic-hexagon.mesh";
+   int ref_levels = 2;
+   int order = 3;
+   bool visualization = 1;
+   double dt = 0.005;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&ref_levels, "-r", "--refine",
+                  "Number of times to refine the mesh uniformly, -1 for auto.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) >= 0.");
+   args.AddOption(&dt, "-dt", "--time-step",
+                  "Time step.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   // 2. Read the mesh from the given mesh file. We can handle triangular,
+   //    quadrilateral, tetrahedral and hexahedral meshes with the same code.
+   Mesh mesh(mesh_file, 1, 1);
+   const int dim = mesh.Dimension();
+
+   // 3. Refine the mesh to increase the resolution. In this example we do
+   //    'ref_levels' of uniform refinement.
+   for (int l = 0; l < ref_levels; l++)
+   {
+      mesh.UniformRefinement();
+   }
+
+   mesh.GetBoundingBox(bb_min, bb_max, max(order, 1));
+
+   // 4. Define a finite element space on the mesh. Here we use discontinuous
+   //    finite elements of the specified order >= 0.
+   DG_FECollection fec(order, dim);
+   FiniteElementSpace fes(&mesh, &fec);
+   cout << "Number of unknowns: " << fes.GetVSize() << endl;
+
+   PartialAssembly dgpa(&fes);
+   Mass mass(&dgpa);
+   MassInverse massinv(&mass);
+   Advection adv_pa(&dgpa, dim);
+
+   // 5. 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);
+   FunctionCoefficient u0(u0_function);
+   u.ProjectCoefficient(u0);
+   {
+      ofstream omesh("ex9.mesh");
+      omesh.precision(8);
+      mesh.Print(omesh);
+      ofstream osol("ex9-init.gf");
+      osol.precision(8);
+      u.Save(osol);
+   }
+
+   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(8);
+         sout << "solution\n" << mesh << u;
+         sout << "pause\n";
+         sout << flush;
+         cout << "GLVis visualization paused."
+              << " Press space (in the GLVis window) to resume it.\n";
+      }
+   }
+
+   // 6. 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).
+   PAEvolution<Advection> adv(&massinv, &adv_pa);
+   ODESolver *ode_solver = new RK4Solver;
+
+   double t = 0.0;
+   double t_final = 10.0;
+   int vis_steps = 5;
+   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;
+         }
+      }
+   }
+
+   // 7. 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(8);
+      u.Save(osol);
+   }
+
+   // 8. Free the used memory.
+   delete ode_solver;
+
+   return 0;
+}
+
+// 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]);
+   }
+   int problem = 0;
+   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;
+            }
+         }
+      }
+      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;
+}

miniapps/conslaw/apps/advection.hpp

@@ -0,0 +1,63 @@
+#ifndef DGPA_ADVECTION
+#define DGPA_ADVECTION
+
+#include "dg.hpp"
+#include "cons_law.hpp"
+
+namespace mfem
+{
+namespace dg
+{
+
+struct DAdv
+{
+   int dim;
+
+   DAdv() : dim(0) { }
+   int NComponents() const
+   {
+      return 1;
+   }
+   void operator()(double *u, double *F) const
+   {
+      for (int d = 0; d < dim; ++d)
+      {
+         F[d] = u[0];
+      }
+   }
+};
+
+struct DAdvUpwinding
+{
+   int dim;
+   DAdvUpwinding() : dim(0) { }
+   int NComponents() const
+   {
+      return 1;
+   }
+   void operator()(double *uL, double *uR, double *n, double *Fhat) const
+   {
+      double bDotN = 0.0;
+      for (int d = 0; d < dim; ++d)
+      {
+         bDotN += 1.0*n[d];
+      }
+      Fhat[0] = 0.0;
+      double u = (bDotN >= 0.0) ? uL[0] : uR[0];
+      for (int d = 0; d < dim; ++d)
+      {
+         Fhat[0] += u*n[d];
+      }
+   }
+};
+
+struct Advection : ConservationLaw<DAdv, DAdvUpwinding>
+{
+   Advection(const PartialAssembly *pa_, int dim)
+      : ConservationLaw(pa_, dim, 1) { }
+};
+
+}
+}
+
+#endif

miniapps/conslaw/apps/burgers.cpp

@@ -0,0 +1,163 @@
+//          MFEM DG FOR CONSERVATION LAWS WITH PARTIAL ASSEMBLY
+
+#include "mfem.hpp"
+#include "dg.hpp"
+#include <fstream>
+#include <iostream>
+#include <random>
+
+#include "burgers.hpp"
+#include "evol.hpp"
+
+using namespace std;
+using namespace mfem;
+using namespace dg;
+
+
+Vector bb_min, bb_max;
+double u0_function(const Vector &x);
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   const char *mesh_file =
+      "../../data/periodic-segment.mesh";
+   int ref_levels = 2;
+   int order = 3;
+   bool visualization = 1;
+   double dt = 0.005;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&ref_levels, "-r", "--refine",
+                  "Number of times to refine the mesh uniformly, -1 for auto.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) >= 0.");
+   args.AddOption(&dt, "-dt", "--time-step",
+                  "Time step.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   // 2. Read the mesh from the given mesh file. We can handle triangular,
+   //    quadrilateral, tetrahedral and hexahedral meshes with the same code.
+   Mesh mesh(mesh_file, 1, 1);
+   const int dim = mesh.Dimension();
+
+   // 3. Refine the mesh to increase the resolution. In this example we do
+   //    'ref_levels' of uniform refinement.
+   for (int l = 0; l < ref_levels; l++)
+   {
+      mesh.UniformRefinement();
+   }
+
+   mesh.GetBoundingBox(bb_min, bb_max, max(order, 1));
+
+   // 4. Define a finite element space on the mesh. Here we use discontinuous
+   //    finite elements of the specified order >= 0.
+   DG_FECollection fec(order, dim);
+   FiniteElementSpace fes(&mesh, &fec);
+   cout << "Number of unknowns: " << fes.GetVSize() << endl;
+
+   PartialAssembly dgpa(&fes);
+   Mass mass(&dgpa);
+   MassInverse massinv(&mass);
+   Burgers burgers(&dgpa, dim, 1);
+
+   // 5. 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);
+   FunctionCoefficient u0(u0_function);
+   u.ProjectCoefficient(u0);
+   {
+      ofstream omesh("ex9.mesh");
+      omesh.precision(8);
+      mesh.Print(omesh);
+      ofstream osol("ex9-init.gf");
+      osol.precision(8);
+      u.Save(osol);
+   }
+
+   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(8);
+         sout << "solution\n" << mesh << u;
+         sout << "pause\n";
+         sout << flush;
+         cout << "GLVis visualization paused."
+              << " Press space (in the GLVis window) to resume it.\n";
+      }
+   }
+
+   // 6. 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).
+   PAEvolution<Burgers> adv(&massinv, &burgers);
+   ODESolver *ode_solver = new RK4Solver;
+
+   double t = 0.0;
+   double t_final = 1.5;
+   int vis_steps = 5;
+   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;
+         }
+      }
+   }
+
+   // 7. 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(8);
+      u.Save(osol);
+   }
+
+   // 8. Free the used memory.
+   delete ode_solver;
+
+   return 0;
+}
+
+// Initial condition
+double u0_function(const Vector &x)
+{
+   return 0.5 + sin(2*M_PI*x(0));
+}

miniapps/conslaw/apps/burgers.hpp

@@ -0,0 +1,62 @@
+#ifndef DGPA_BURGERS
+#define DGPA_BURGERS
+
+#include "dg.hpp"
+#include "cons_law.hpp"
+
+namespace mfem
+{
+namespace dg
+{
+
+struct DBurgers
+{
+   int dim;
+
+   DBurgers() : dim(0) { }
+   int NComponents() const
+   {
+      return 1;
+   }
+   void operator()(double *u, double *F) const
+   {
+      for (int d = 0; d < dim; ++d)
+      {
+         F[d] = 0.5*u[0]*u[0];
+      }
+   }
+};
+
+struct DBurgersRiemann
+{
+   int dim;
+   DBurgersRiemann() : dim(0) { }
+   int NComponents() const
+   {
+      return 1;
+   }
+   void operator()(double *uL, double *uR, double *n, double *Fhat) const
+   {
+      double cL = 0.5*uL[0]*n[0];
+      double cR = 0.5*uR[0]*n[0];
+      double avg = 0.5*(cL + cR);
+
+      if (cL > cR)
+      {
+         Fhat[0] = (avg > 0) ? cL*uL[0] : cR*uR[0];
+      }
+      else
+      {
+         if (cL > 0) { Fhat[0] = cL*uL[0]; }
+         else if (cL*cR < 0) { Fhat[0] = 0; }
+         else { Fhat[0] = cR*uR[0]; }
+      }
+   }
+};
+
+using Burgers = ConservationLaw<DBurgers, DBurgersRiemann>;
+
+}
+}
+
+#endif

miniapps/conslaw/apps/cons_law.hpp

@@ -0,0 +1,40 @@
+#ifndef DGPA_CONS_LAW
+#define DGPA_CONS_LAW
+
+#include "dg.hpp"
+
+namespace mfem
+{
+namespace dg
+{
+
+template <typename F, typename Fhat>
+class ConservationLaw : Operator
+{
+   const PartialAssembly *pa;
+   GtDB<F> vol;
+   BtDB_face<Fhat> face;
+   int nc;
+public:
+   ConservationLaw(const PartialAssembly *pa_, int dim, int nc_)
+      : pa(pa_), vol(pa), face(pa), nc(nc_)
+   {
+      vol.d.dim = dim;
+      face.d.dim = dim;
+   }
+   void Mult(const Vector &x, Vector &y) const
+   {
+      y = 0.0;
+      vol.Mult(x, y);
+      face.Mult(x, y);
+   }
+   int Size() const
+   {
+      return pa->GetFES()->GetNDofs()*nc;
+   }
+};
+
+}
+}
+
+#endif

miniapps/conslaw/apps/euler.cpp

@@ -0,0 +1,260 @@
+//          MFEM DG FOR CONSERVATION LAWS WITH PARTIAL ASSEMBLY
+
+#include "mfem.hpp"
+#include "dg.hpp"
+#include <fstream>
+#include <iostream>
+#include <random>
+
+#include "euler.hpp"
+#include "evol.hpp"
+
+using namespace std;
+using namespace mfem;
+using namespace dg;
+
+int problem = 1;
+int precision = 8;
+
+void InitialCondition(const Vector &x, Vector &y)
+{
+   MFEM_ASSERT(x.Size() == 2, "");
+
+   double radius = 0, Minf = 0, beta = 0;
+   if (problem == 1)
+   {
+      // "Fast vortex"
+      radius = 0.2;
+      Minf = 0.5;
+      beta = 1. / 5.;
+   }
+   else if (problem == 2)
+   {
+      // "Slow vortex"
+      radius = 0.2;
+      Minf = 0.05;
+      beta = 1. / 50.;
+   }
+   else
+   {
+      mfem_error("Cannot recognize problem."
+                 "Options are: 1 - fast vortex, 2 - slow vortex");
+   }
+
+   const double specific_heat_ratio = 1.4;
+   const double gas_constant = 1.0;
+
+   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;
+}
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   const char *mesh_file =
+      "../../data/periodic-square.mesh";
+   int ref_levels = 1;
+   int order = 2;
+   bool visualization = 1;
+   double dt = 0.005;
+   double t_final = 2.0;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&ref_levels, "-r", "--refine",
+                  "Number of times to refine the mesh uniformly, -1 for auto.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) >= 0.");
+   args.AddOption(&dt, "-dt", "--time-step",
+                  "Time step.");
+   args.AddOption(&t_final, "-tf", "--t-final",
+                  "Final time; start time is 0.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   // 2. Read the mesh from the given mesh file. We can handle triangular,
+   //    quadrilateral, tetrahedral and hexahedral meshes with the same code.
+   //    NURBS meshes are projected to second order meshes.
+   Mesh mesh(mesh_file, 1, 1);
+   const int dim = mesh.Dimension();
+   const int num_equation = dim + 2;
+
+   // 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();
+   }
+
+   // 4. Define a finite element space on the mesh. Here we use discontinuous
+   //    finite elements of the specified order >= 0.
+   DG_FECollection fec(order, dim);
+   FiniteElementSpace scalar_fes(&mesh, &fec);
+   FiniteElementSpace momentum_fes(&mesh, &fec, dim, Ordering::byNODES);
+   FiniteElementSpace fes(&mesh, &fec, num_equation, Ordering::byNODES);
+   cout << "Number of unknowns: " << fes.GetVSize() << endl;
+
+   PartialAssembly dgpa(&fes);
+   Mass mass(&dgpa);
+   MassInverse massinv(&mass);
+
+   Array<int> offsets(num_equation + 1);
+   for (int k = 0; k <= num_equation; k++) { offsets[k] = k * fes.GetNDofs(); }
+   BlockVector u_block(offsets);
+
+   VectorFunctionCoefficient u0(num_equation, InitialCondition);
+   GridFunction u(&fes, u_block.GetData());
+   GridFunction rho(&scalar_fes, u_block.GetData());
+   GridFunction rho_u(&momentum_fes, u_block.GetData() + offsets[1]);
+   u.ProjectCoefficient(u0);
+
+   Euler euler(&dgpa, dim);
+
+   // Output the initial solution.
+   {
+      ofstream mesh_ofs("vortex.mesh");
+      mesh_ofs.precision(precision);
+      mesh_ofs << mesh;
+
+      for (int k = 0; k < num_equation; k++)
+      {
+         GridFunction uk(&scalar_fes, u_block.GetBlock(k));
+         ostringstream sol_name;
+         sol_name << "vortex-" << k << "-init.gf";
+         ofstream sol_ofs(sol_name.str().c_str());
+         sol_ofs.precision(precision);
+         sol_ofs << uk;
+      }
+   }
+
+   // 5. 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).
+   PAEvolution<Euler> euler_evol(&massinv, &euler);
+   ODESolver *ode_solver = new RK4Solver;
+
+   // Visualize the density
+   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 << rho_u;
+         sout << "pause\n";
+         sout << flush;
+         cout << "GLVis visualization paused."
+              << " Press space (in the GLVis window) to resume it.\n";
+      }
+   }
+
+   // Start the timer.
+   tic_toc.Clear();
+   tic_toc.Start();
+
+   int vis_steps = 50;
+   double t = 0.0;
+   euler_evol.SetTime(t);
+   ode_solver->Init(euler_evol);
+
+   // Integrate in time.
+   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 << rho_u << flush;
+         }
+      }
+   }
+
+   tic_toc.Stop();
+   cout << " done, " << tic_toc.RealTime() << "s." << endl;
+
+   // 6. 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++)
+   {
+      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;
+   }
+
+   // 7. Compute the L2 solution error summed for all components.
+   if (t_final == 2.0)
+   {
+      const double error = u.ComputeLpError(2, u0);
+      cout << "Solution error: " << error << endl;
+   }
+
+   // 8. Free the used memory.
+   delete ode_solver;
+
+   return 0;
+}

miniapps/conslaw/apps/euler.hpp

@@ -0,0 +1,47 @@
+#ifndef DG_EULER
+#define DG_EULER
+
+#include "cons_law.hpp"
+#include "ns_autogen.hpp"
+
+namespace mfem
+{
+namespace dg
+{
+
+struct EulerFlux
+{
+   int dim = 0;
+   int NComponents() const
+   {
+      return dim + 2;
+   }
+   void operator()(double *u, double *F) const
+   {
+      eulerF(dim, u, F);
+   }
+};
+
+struct EulerNumericalFlux
+{
+   int dim = 0;
+   int NComponents() const
+   {
+      return dim + 2;
+   }
+   void operator()(double *uL, double *uR, double *n, double *Fhat) const
+   {
+      eulerFhat(dim, uR, uL, n, Fhat);
+   }
+};
+
+struct Euler : ConservationLaw<EulerFlux, EulerNumericalFlux>
+{
+   Euler(const PartialAssembly *pa_, int dim)
+      : ConservationLaw(pa_, dim, dim + 2) { }
+};
+
+}
+}
+
+#endif

miniapps/conslaw/apps/evol.hpp

@@ -0,0 +1,35 @@
+#ifndef DGPA_EVOL
+#define DGPA_EVOL
+
+#include "dg.hpp"
+
+namespace mfem
+{
+namespace dg
+{
+
+template <typename Oper>
+class PAEvolution : public TimeDependentOperator
+{
+private:
+   MassInverse *Minv;
+   Oper *oper;
+   mutable Vector z;
+public:
+   PAEvolution(MassInverse *Minv_, Oper *oper_)
+      : TimeDependentOperator(oper_->Size()), Minv(Minv_), oper(oper_) { }
+
+   virtual void Mult(const Vector &x, Vector &y) const
+   {
+      z.SetSize(x.Size());
+      oper->Mult(x, z);
+      Minv->Mult(z, y);
+   }
+
+   virtual ~PAEvolution() { }
+};
+
+}
+}
+
+#endif

miniapps/conslaw/apps/mass.hpp

@@ -0,0 +1,39 @@
+#ifndef DG_MASS_OPER
+#define DG_MASS_OPER
+
+#include "dg.hpp"
+
+namespace mfem
+{
+namespace dg
+{
+
+struct DId
+{
+   int NComponents() const
+   {
+      return 1;
+   }
+   void operator()(double *u, double *F) const
+   {
+      F[0] = u[0];
+   }
+};
+
+class DGMassPA : Operator
+{
+   const PartialAssembly *pa;
+   BtDB<DId> oper;
+public:
+   DGMassPA(const PartialAssembly *pa_) : pa(pa_), oper(pa) { }
+   void Mult(const Vector &x, Vector &y) const
+   {
+      y = 0.0;
+      oper.Mult(x, y);
+   }
+};
+
+}
+}
+
+#endif

miniapps/conslaw/apps/ns_autogen.hpp

@@ -0,0 +1,567 @@
+#include <math.h>
+
+/******************************************************************************/
+
+void roe1d (
+   double UR[3],
+   double UL[3],
+   double n[1],
+   double F[3])
+{
+   double t76;
+   double t64;
+   double t65;
+   double t46;
+   double t21;
+   double t22;
+   double t23;
+   double t5;
+   double t100;
+   double t29;
+   double t8;
+   double t1;
+   double t2;
+   double t30;
+   double t9;
+   double t10;
+   double t24;
+   double t25;
+   double t26;
+   double t79;
+   double t67;
+   double t68;
+   double t70;
+   double t47;
+   double t48;
+   double t49;
+   double t12;
+   double t14;
+   double t107;
+   double t88;
+   double t89;
+   double t19;
+   double t43;
+   double t86;
+   double t62;
+   double t81;
+   double t53;
+   double t35;
+   double t36;
+   double t56;
+   double t33;
+   double t34;
+   double t57;
+   double t59;
+   double t61;
+   double t73;
+   t1 = UR[1];
+   t2 = n[0];
+   t5 = UL[1];
+   t8 = UR[0];
+   t9 = UL[0];
+   t10 = 0.1e1 / t9;
+   t12 = sqrt(t8 * t10);
+   t14 = 0.1e1 / t8;
+   t19 = 0.100e1 * t12 * t1 * t14 + 0.10e1 * t5 * t10;
+   t21 = 0.10e1 * t12 + 0.10e1;
+   t22 = 0.1e1 / t21;
+   t23 = t19 * t22;
+   t24 = t23 * t2;
+   t25 = fabs(t24);
+   t26 = t8 - t9;
+   t29 = 0.10e1 * t24;
+   t30 = UR[2];
+   t33 = 0.4e0 * t30;
+   t34 = t1 * t1;
+   t35 = t14 * t34;
+   t36 = 0.2000e0 * t35;
+   t43 = UL[2];
+   t46 = 0.4e0 * t43;
+   t47 = t5 * t5;
+   t48 = t10 * t47;
+   t49 = 0.2000e0 * t48;
+   t53 = 0.10e1 * t12 * (0.10e1 * t30 * t14 + 0.10e1 * (t33 - t36) * t14) + 0.10e1
+         * t43 * t10 + 0.10e1 * (t46 - t49) * t10;
+   t56 = t19 * t19;
+   t57 = t21 * t21;
+   t59 = t56 / t57;
+   t61 = 0.40e0 * t53 * t22 - 0.2000e0 * t59;
+   t62 = sqrt(t61);
+   t64 = fabs(t29 + t62);
+   t65 = 0.5e0 * t64;
+   t67 = fabs(-t29 + t62);
+   t68 = 0.5e0 * t67;
+   t70 = t65 + t68 - 0.10e1 * t25;
+   t73 = t1 - t5;
+   t76 = 0.2000e0 * t59 * t26 - 0.40e0 * t23 * t73 + t33 - t46;
+   t79 = t70 * t76 / t61;
+   t81 = t65 - t68;
+   t86 = -0.10e1 * t23 * t2 * t26 + t73 * t2;
+   t88 = 0.1e1 / t62;
+   t89 = t81 * t86 * t88;
+   F[0] = 0.50e0 * t1 * t2 + 0.50e0 * t5 * t2 - 0.50e0 * t25 * t26 - 0.5e0 * t79 -
+          0.5e0 * t89;
+   t100 = t79 + t89;
+   t107 = t81 * t76 * t88 + t70 * t86;
+   F[1] = 0.50e0 * t35 * t2 + 0.50e0 * t48 * t2 + 0.5e0 * t2 *
+          (t33 - t36 + t46 - t49) - 0.50e0 * t25 * t73 - 0.50e0 * t100 * t19 * t22 - 0.5e0
+          * t107 * t2;
+   F[2] = 0.50e0 * (0.14e1 * t30 - t36) * t1 * t14 * t2 + 0.50e0 *
+          (0.14e1 * t43 - t49) * t5 * t10 * t2 - 0.50e0 * t25 * (t30 - t43) - 0.50e0 *
+          t100 * t53 * t22 - 0.50e0 * t107 * t19 * t22 * t2;
+   return;
+}
+
+/******************************************************************************/
+
+void roe2d (
+   double UR[4],
+   double UL[4],
+   double n[2],
+   double F[4])
+{
+   double t17;
+   double t46;
+   double t155;
+   double t23;
+   double t59;
+   double t35;
+   double t76;
+   double t47;
+   double t79;
+   double t49;
+   double t50;
+   double t51;
+   double t80;
+   double t2;
+   double t13;
+   double t81;
+   double t82;
+   double t22;
+   double t98;
+   double t99;
+   double t8;
+   double t85;
+   double t52;
+   double t55;
+   double t58;
+   double t37;
+   double t38;
+   double t39;
+   double t16;
+   double t69;
+   double t100;
+   double t90;
+   double t106;
+   double t105;
+   double t103;
+   double t102;
+   double t109;
+   double t108;
+   double t3;
+   double t5;
+   double t60;
+   double t61;
+   double t112;
+   double t111;
+   double t113;
+   double t94;
+   double t119;
+   double t122;
+   double t18;
+   double t19;
+   double t125;
+   double t41;
+   double t30;
+   double t97;
+   double t128;
+   double t130;
+   double t26;
+   double t64;
+   double t10;
+   double t134;
+   double t1;
+   double t137;
+   double t136;
+   double t143;
+   double t148;
+   t1 = UR[0];
+   t2 = UR[1];
+   t3 = 0.1e1 / t1;
+   t5 = n[0];
+   t8 = UR[2];
+   t10 = n[1];
+   t13 = 0.10e1 * t2 * t3 * t5 + 0.10e1 * t8 * t3 * t10;
+   t16 = UL[0];
+   t17 = UL[1];
+   t18 = 0.1e1 / t16;
+   t19 = t17 * t18;
+   t22 = UL[2];
+   t23 = t22 * t18;
+   t26 = 0.10e1 * t19 * t5 + 0.10e1 * t23 * t10;
+   t30 = sqrt(t1 * t18);
+   t35 = 0.100e1 * t30 * t2 * t3 + 0.10e1 * t19;
+   t37 = 0.10e1 * t30 + 0.10e1;
+   t38 = 0.1e1 / t37;
+   t39 = t35 * t38;
+   t41 = 0.10e1 * t39 * t5;
+   t46 = 0.100e1 * t30 * t8 * t3 + 0.10e1 * t23;
+   t47 = t46 * t38;
+   t49 = 0.10e1 * t47 * t10;
+   t50 = t41 + t49;
+   t51 = fabs(t50);
+   t52 = t1 - t16;
+   t55 = UR[3];
+   t58 = 0.4e0 * t55;
+   t59 = t2 * t2;
+   t60 = t1 * t1;
+   t61 = 0.1e1 / t60;
+   t64 = t8 * t8;
+   t69 = 0.20e0 * t1 * (0.100e1 * t59 * t61 + 0.100e1 * t64 * t61);
+   t76 = UL[3];
+   t79 = 0.4e0 * t76;
+   t80 = t17 * t17;
+   t81 = t16 * t16;
+   t82 = 0.1e1 / t81;
+   t85 = t22 * t22;
+   t90 = 0.20e0 * t16 * (0.100e1 * t80 * t82 + 0.100e1 * t85 * t82);
+   t94 = 0.10e1 * t30 * (0.10e1 * t55 * t3 + 0.10e1 * (t58 - t69) * t3) + 0.10e1 *
+         t76 * t18 + 0.10e1 * (t79 - t90) * t18;
+   t97 = t35 * t35;
+   t98 = t37 * t37;
+   t99 = 0.1e1 / t98;
+   t100 = t97 * t99;
+   t102 = t46 * t46;
+   t103 = t102 * t99;
+   t105 = 0.40e0 * t94 * t38 - 0.2000e0 * t100 - 0.2000e0 * t103;
+   t106 = sqrt(t105);
+   t108 = fabs(t41 + t49 + t106);
+   t109 = 0.5e0 * t108;
+   t111 = fabs(-t41 - t49 + t106);
+   t112 = 0.5e0 * t111;
+   t113 = t109 + t112 - t51;
+   t119 = t2 - t17;
+   t122 = t8 - t22;
+   t125 = 0.4e0 * (0.500e0 * t100 + 0.500e0 * t103) * t52 - 0.40e0 * t39 * t119 -
+          0.40e0 * t47 * t122 + t58 - t79;
+   t128 = t113 * t125 / t105;
+   t130 = t109 - t112;
+   t134 = -t50 * t52 + t119 * t5 + t122 * t10;
+   t136 = 0.1e1 / t106;
+   t137 = t130 * t134 * t136;
+   F[0] = 0.5e0 * t1 * t13 + 0.5e0 * t16 * t26 - 0.5e0 * t51 * t52 - 0.5e0 * t128 -
+          0.5e0 * t137;
+   t143 = t58 - t69 + t79 - t90;
+   t148 = t128 + t137;
+   t155 = t130 * t125 * t136 + t113 * t134;
+   F[1] = 0.5e0 * t2 * t13 + 0.5e0 * t17 * t26 + 0.5e0 * t5 * t143 - 0.5e0 * t51 *
+          t119 - 0.50e0 * t148 * t35 * t38 - 0.5e0 * t155 * t5;
+   F[2] = 0.5e0 * t8 * t13 + 0.5e0 * t22 * t26 + 0.5e0 * t10 * t143 - 0.5e0 * t51 *
+          t122 - 0.50e0 * t148 * t46 * t38 - 0.5e0 * t155 * t10;
+   F[3] = 0.5e0 * (0.14e1 * t55 - t69) * t13 + 0.5e0 * (0.14e1 * t76 - t90) * t26 -
+          0.5e0 * t51 * (t55 - t76) - 0.50e0 * t148 * t94 * t38 - 0.5e0 * t155 * t50;
+   return;
+}
+
+/******************************************************************************/
+
+void roe3d (
+   double UR[5],
+   double UL[5],
+   double n[3],
+   double F[5])
+{
+   double t105;
+   double t138;
+   double t8;
+   double t10;
+   double t44;
+   double t27;
+   double t28;
+   double t35;
+   double t186;
+   double t76;
+   double t55;
+   double t56;
+   double t58;
+   double t23;
+   double t24;
+   double t18;
+   double t125;
+   double t126;
+   double t128;
+   double t129;
+   double t131;
+   double t132;
+   double t134;
+   double t135;
+   double t137;
+   double t77;
+   double t78;
+   double t81;
+   double t48;
+   double t108;
+   double t113;
+   double t117;
+   double t120;
+   double t121;
+   double t122;
+   double t123;
+   double t149;
+   double t46;
+   double t47;
+   double t50;
+   double t63;
+   double t64;
+   double t2;
+   double t3;
+   double t152;
+   double t155;
+   double t158;
+   double t160;
+   double t96;
+   double t66;
+   double t67;
+   double t68;
+   double t165;
+   double t167;
+   double t168;
+   double t1;
+   double t69;
+   double t174;
+   double t13;
+   double t5;
+   double t89;
+   double t139;
+   double t179;
+   double t31;
+   double t32;
+   double t21;
+   double t22;
+   double t84;
+   double t15;
+   double t39;
+   double t75;
+   double t99;
+   double t72;
+   double t146;
+   double t100;
+   double t101;
+   double t102;
+   t1 = UR[0];
+   t2 = UR[1];
+   t3 = 0.1e1 / t1;
+   t5 = n[0];
+   t8 = UR[2];
+   t10 = n[1];
+   t13 = UR[3];
+   t15 = n[2];
+   t18 = 0.10e1 * t2 * t3 * t5 + 0.10e1 * t8 * t3 * t10 + 0.10e1 * t13 * t3 * t15;
+   t21 = UL[0];
+   t22 = UL[1];
+   t23 = 0.1e1 / t21;
+   t24 = t22 * t23;
+   t27 = UL[2];
+   t28 = t27 * t23;
+   t31 = UL[3];
+   t32 = t31 * t23;
+   t35 = 0.10e1 * t24 * t5 + 0.10e1 * t10 * t28 + 0.10e1 * t32 * t15;
+   t39 = sqrt(t1 * t23);
+   t44 = 0.100e1 * t39 * t2 * t3 + 0.10e1 * t24;
+   t46 = 0.10e1 * t39 + 0.10e1;
+   t47 = 0.1e1 / t46;
+   t48 = t44 * t47;
+   t50 = 0.10e1 * t48 * t5;
+   t55 = 0.100e1 * t39 * t8 * t3 + 0.10e1 * t28;
+   t56 = t55 * t47;
+   t58 = 0.10e1 * t56 * t10;
+   t63 = 0.100e1 * t39 * t13 * t3 + 0.10e1 * t32;
+   t64 = t63 * t47;
+   t66 = 0.10e1 * t64 * t15;
+   t67 = t50 + t58 + t66;
+   t68 = fabs(t67);
+   t69 = t1 - t21;
+   t72 = UR[4];
+   t75 = 0.4e0 * t72;
+   t76 = t2 * t2;
+   t77 = t1 * t1;
+   t78 = 0.1e1 / t77;
+   t81 = t8 * t8;
+   t84 = t13 * t13;
+   t89 = 0.20e0 * t1 * (0.100e1 * t76 * t78 + 0.100e1 * t81 * t78 + 0.100e1 * t84 *
+                        t78);
+   t96 = UL[4];
+   t99 = 0.4e0 * t96;
+   t100 = t22 * t22;
+   t101 = t21 * t21;
+   t102 = 0.1e1 / t101;
+   t105 = t27 * t27;
+   t108 = t31 * t31;
+   t113 = 0.20e0 * t21 * (0.100e1 * t100 * t102 + 0.100e1 * t105 * t102 + 0.100e1 *
+                          t108 * t102);
+   t117 = 0.10e1 * t39 * (0.10e1 * t72 * t3 + 0.10e1 * (t75 - t89) * t3) + 0.10e1 *
+          t96 * t23 + 0.10e1 * (t99 - t113) * t23;
+   t120 = t44 * t44;
+   t121 = t46 * t46;
+   t122 = 0.1e1 / t121;
+   t123 = t120 * t122;
+   t125 = t55 * t55;
+   t126 = t125 * t122;
+   t128 = t63 * t63;
+   t129 = t128 * t122;
+   t131 = 0.40e0 * t117 * t47 - 0.2000e0 * t123 - 0.2000e0 * t126 - 0.2000e0 *
+          t129;
+   t132 = sqrt(t131);
+   t134 = fabs(t50 + t58 + t66 + t132);
+   t135 = 0.5e0 * t134;
+   t137 = fabs(-t50 - t58 - t66 + t132);
+   t138 = 0.5e0 * t137;
+   t139 = t135 + t138 - t68;
+   t146 = t2 - t22;
+   t149 = t8 - t27;
+   t152 = t13 - t31;
+   t155 = 0.4e0 * (0.500e0 * t123 + 0.500e0 * t126 + 0.500e0 * t129) * t69 - 0.40e0
+          * t48 * t146 - 0.40e0 * t56 * t149 - 0.40e0 * t64 * t152 + t75 - t99;
+   t158 = t139 * t155 / t131;
+   t160 = t135 - t138;
+   t165 = -t67 * t69 + t146 * t5 + t149 * t10 + t152 * t15;
+   t167 = 0.1e1 / t132;
+   t168 = t160 * t165 * t167;
+   F[0] = 0.5e0 * t1 * t18 + 0.5e0 * t21 * t35 - 0.5e0 * t68 * t69 - 0.5e0 * t158 -
+          0.5e0 * t168;
+   t174 = t75 - t89 + t99 - t113;
+   t179 = t158 + t168;
+   t186 = t160 * t155 * t167 + t139 * t165;
+   F[1] = 0.5e0 * t2 * t18 + 0.5e0 * t22 * t35 + 0.5e0 * t5 * t174 - 0.5e0 * t68 *
+          t146 - 0.50e0 * t179 * t44 * t47 - 0.5e0 * t186 * t5;
+   F[2] = 0.5e0 * t8 * t18 + 0.5e0 * t27 * t35 + 0.5e0 * t10 * t174 - 0.5e0 * t68 *
+          t149 - 0.50e0 * t179 * t55 * t47 - 0.5e0 * t186 * t10;
+   F[3] = 0.5e0 * t13 * t18 + 0.5e0 * t31 * t35 + 0.5e0 * t15 * t174 - 0.5e0 * t68
+          * t152 - 0.50e0 * t179 * t63 * t47 - 0.5e0 * t186 * t15;
+   F[4] = 0.5e0 * (0.14e1 * t72 - t89) * t18 + 0.5e0 * (0.14e1 * t96 - t113) * t35
+          - 0.5e0 * t68 * (t72 - t96) - 0.50e0 * t179 * t117 * t47 - 0.5e0 * t186 * t67;
+   return;
+}
+
+/******************************************************************************/
+
+void Fi1d (double U[3], double Fi[3])
+{
+   double t4;
+   double t1;
+   double t3;
+   double t6;
+   Fi[0] = U[1];
+   t1 = pow(Fi[0], 0.2e1);
+   t3 = 0.1e1 / U[0];
+   t4 = t1 * t3;
+   t6 = U[2];
+   Fi[1] = 0.4e1 / 0.5e1 * t4 + 0.2e1 / 0.5e1 * t6;
+   Fi[2] = Fi[0] * (0.7e1 / 0.5e1 * t6 - t4 / 0.5e1) * t3;
+   return;
+}
+
+/******************************************************************************/
+
+void Fi2d (double U[4], double Fi[8])
+{
+   double t1;
+   double t3;
+   double t5;
+   double t6;
+   double t7;
+   double t8;
+   double t11;
+   double t14;
+   Fi[0] = U[1];
+   t1 = pow(Fi[0], 0.2e1);
+   t3 = 0.1e1 / U[0];
+   t5 = U[3];
+   t6 = 0.2e1 / 0.5e1 * t5;
+   t7 = U[2];
+   t8 = t7 * t7;
+   t11 = (t1 + t8) * t3 / 0.5e1;
+   Fi[1] = t1 * t3 + t6 - t11;
+   Fi[2] = Fi[0] * t7 * t3;
+   t14 = 0.7e1 / 0.5e1 * t5 - t11;
+   Fi[3] = Fi[0] * t14 * t3;
+   Fi[4] = t7;
+   Fi[5] = Fi[2];
+   Fi[6] = t8 * t3 + t6 - t11;
+   Fi[7] = Fi[4] * t14 * t3;
+   return;
+}
+
+/******************************************************************************/
+
+void Fi3d (double U[5], double Fi[15])
+{
+   double t17;
+   double t10;
+   double t9;
+   double t13;
+   double t3;
+   double t5;
+   double t6;
+   double t1;
+   double t7;
+   double t8;
+   Fi[0] = U[1];
+   t1 = pow(Fi[0], 0.2e1);
+   t3 = 0.1e1 / U[0];
+   t5 = U[4];
+   t6 = 0.2e1 / 0.5e1 * t5;
+   t7 = U[2];
+   t8 = t7 * t7;
+   t9 = U[3];
+   t10 = t9 * t9;
+   t13 = (t1 + t8 + t10) * t3 / 0.5e1;
+   Fi[1] = t1 * t3 + t6 - t13;
+   Fi[2] = Fi[0] * t7 * t3;
+   Fi[3] = Fi[0] * t9 * t3;
+   t17 = 0.7e1 / 0.5e1 * t5 - t13;
+   Fi[4] = Fi[0] * t17 * t3;
+   Fi[5] = t7;
+   Fi[6] = Fi[2];
+   Fi[7] = t8 * t3 + t6 - t13;
+   Fi[8] = Fi[5] * t9 * t3;
+   Fi[9] = Fi[5] * t17 * t3;
+   Fi[10] = t9;
+   Fi[11] = Fi[3];
+   Fi[12] = Fi[8];
+   Fi[13] = t10 * t3 + t6 - t13;
+   Fi[14] = Fi[10] * t17 * t3;
+   return;
+}
+
+/******************************************************************************/
+void eulerF(int dim, double *U, double *F)
+{
+   if (dim == 1)
+   {
+      Fi1d(U,F);
+   }
+   else if (dim == 2)
+   {
+      Fi2d(U,F);
+   }
+   else if (dim == 3)
+   {
+      Fi3d(U,F);
+   }
+}
+void eulerFhat(int dim, double *UR, double *UL, double *n, double *F)
+{
+   if (dim == 1)
+   {
+      roe1d(UR,UL,n,F);
+   }
+   else if (dim == 2)
+   {
+      roe2d(UR,UL,n,F);
+   }
+   else if (dim == 3)
+   {
+      roe3d(UR,UL,n,F);
+   }
+}

miniapps/conslaw/apps/tests.cpp

@@ -0,0 +1,122 @@
+#include "mfem.hpp"
+#include "advection.hpp"
+#include "mass.hpp"
+#include <fstream>
+#include <iostream>
+#include <random>
+
+using namespace std;
+using namespace mfem;
+using namespace dg;
+
+void velocity_function(const Vector &x, Vector &v)
+{
+   int dim = x.Size();
+   for (int d = 0; d < dim; ++d)
+   {
+      v(d) = 1.0;
+   }
+}
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   const char *mesh_file =
+      "../../data/periodic-hexagon.mesh";
+   int ref_levels = -1;
+   int order = 1;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&ref_levels, "-r", "--refine",
+                  "Number of times to refine the mesh uniformly, -1 for auto.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) >= 0.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   // 2. Read the mesh from the given mesh file. We can handle triangular,
+   //    quadrilateral, tetrahedral and hexahedral meshes with the same code.
+
+   Mesh mesh(mesh_file, 1, 1);
+   const int dim = mesh.Dimension();
+
+   // 3. Refine the mesh to increase the resolution. In this example we do
+   //    'ref_levels' of uniform refinement. By default, or if ref_levels < 0,
+   //    we choose it to be the largest number that gives a final mesh with no
+   //    more than 50,000 elements.
+   {
+      if (ref_levels < 0)
+      {
+         ref_levels = (int)floor(log(50000./mesh.GetNE())/log(2.)/dim);
+      }
+      for (int l = 0; l < ref_levels; l++)
+      {
+         mesh.UniformRefinement();
+      }
+   }
+
+   // 4. Define a finite element space on the mesh. Here we use discontinuous
+   //    finite elements of the specified order >= 0.
+   DG_FECollection fec(order, dim);
+   FiniteElementSpace fes(&mesh, &fec);
+   cout << "Number of unknowns: " << fes.GetVSize() << endl;
+
+   PartialAssembly dgpa(&fes);
+   Mass mass(&dgpa);
+   MassInverse massinv(&mass);
+
+   // Test mass matrix routines
+   GridFunction u(&fes), u2(&fes), Mu(&fes), Mu2(&fes);
+   // Generate random grid function
+   mt19937 re(20);
+   uniform_real_distribution<double> unif(0.0,1.0);
+   for (int i = 0; i < u.Size(); ++i) { u[i] = unif(re); }
+
+   mass.Mult(u, Mu);
+   massinv.Mult(Mu, u2);
+   u2 -= u;
+   cout << "Difference u and M^{-1}M u = " << u2.Normlinf() << '\n';
+
+   DGMassPA mass_pa(&dgpa);
+   mass_pa.Mult(u, Mu2);
+   Mu2 -= Mu;
+   cout << "Difference M u and M_{PA} u = " << Mu2.Normlinf() << '\n';
+
+   // Test advection integrators
+   Advection adv_pa(&dgpa, dim);
+   GridFunction Au(&fes), Au2(&fes);
+
+   // Compare with standard MFEM integrator
+   tic();
+   BilinearForm k(&fes);
+   VectorFunctionCoefficient velocity(dim, velocity_function);
+   k.AddDomainIntegrator(new ConvectionIntegrator(velocity, -1.0));
+   k.AddInteriorFaceIntegrator(
+      new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
+   k.AddBdrFaceIntegrator(
+      new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
+   k.Assemble();
+   k.Finalize();
+   toc();
+   std::cout << "Standard assembly:    " << tic_toc.RealTime() << " s\n";
+   tic();
+   k.Mult(u, Au);
+   toc();
+   std::cout << "Standard application: " << tic_toc.RealTime() << " s\n";
+   // Compute using partial assembly
+   tic();
+   adv_pa.Mult(u, Au2);
+   toc();
+   std::cout << "PA application:       " << tic_toc.RealTime() << " s\n";
+   Au2 -= Au;
+   cout << "Difference A u and A_{PA} u = " << Au2.Normlinf() << '\n';
+
+   return 0;
+}

miniapps/conslaw/include/dg.hpp

@@ -0,0 +1,10 @@
+// Convenience header for DG partial assembly includes
+
+#ifndef DG_HPP
+#define DG_HPP
+
+#include "dg_pa.hpp"
+#include "dg_mass.hpp"
+#include "dg_opers.hpp"
+
+#endif

miniapps/conslaw/include/dg_mass.hpp

@@ -0,0 +1,39 @@
+#ifndef DGPA_MASS
+#define DGPA_MASS
+
+#include "dg_pa.hpp"
+
+namespace mfem
+{
+namespace dg
+{
+
+class Mass : public Operator
+{
+   friend class MassInverse;
+
+   const PartialAssembly *pa;
+   const FiniteElementSpace *fes;
+   DenseTensor M;
+public:
+   Mass(const PartialAssembly *pa_);
+   void Mult(const Vector &x, Vector &y) const;
+   const PartialAssembly *GetPA() const;
+};
+
+class MassInverse : public Operator
+{
+   const Mass *mass;
+   const FiniteElementSpace *fes;
+   DenseTensor Minv;
+public:
+   MassInverse(const Mass *mass_);
+   virtual void Mult(const Vector &x, Vector &y) const;
+};
+
+
+
+} // namespace dg
+} // namespace mfem
+
+#endif

miniapps/conslaw/include/dg_metric.hpp

@@ -0,0 +1,36 @@
+#ifndef DGPA_METRIC
+#define DGPA_METRIC
+
+#include "mfem.hpp"
+
+namespace mfem
+{
+namespace dg
+{
+
+class MetricTerms
+{
+   int nel, nquad, dim, nfaces, nquad_face;
+   Vector detJ;
+   Vector w, wface;
+   DenseTensor Jinv;
+   mutable DenseMatrix nvec;
+   mutable Vector n;
+public:
+   MetricTerms();
+   void Precompute(const FiniteElementSpace *fes,
+                   const IntegrationRule *ir,
+                   const IntegrationRule *ir_face);
+   // Interface to access metric terms at quadrature points:
+   double JacobianDeterminant(int elid, int iq) const;
+   double Weight(int elid, int iq) const;
+   const DenseMatrix& InverseJacobian(int elid, int iq) const;
+   const Vector& Normal(int fid, int iq) const;
+   double FaceWeight(int fid, int iq) const;
+};
+
+
+} // namespace dg
+} // namespace mfem
+
+#endif

miniapps/conslaw/include/dg_opers.hpp

@@ -0,0 +1,210 @@
+#ifndef DGPA_OPERS
+#define DGPA_OPERS
+
+#include "dg_pa.hpp"
+
+namespace mfem
+{
+namespace dg
+{
+
+template <typename D>
+class BtDB : public Operator
+{
+   const PartialAssembly *pa;
+   const FiniteElementSpace *fes;
+public:
+   D d;
+   BtDB(const PartialAssembly *pa_) : pa(pa_), fes(pa->GetFES()) { }
+   void Mult(const Vector &x, Vector &y) const
+   {
+      const int nc = d.NComponents();
+      Array<int> vdofs;
+      DenseMatrix xquad, yquad;
+      DenseMatrix xel, yel;
+      Vector xpt(nc), ypt(nc);
+      for (int i = 0; i < fes->GetNE(); i++)
+      {
+         const DenseMatrix &B = pa->BasisEval(i);
+         const int nquad = B.Height();
+         const int ndof = B.Width();
+         xquad.SetSize(nquad, nc);
+         yquad.SetSize(nquad, nc);
+         xel.SetSize(ndof, nc);
+         yel.SetSize(ndof, nc);
+         fes->GetElementVDofs(i, vdofs);
+         x.GetSubVector(vdofs, xel.Data());
+         mfem::Mult(B, xel, xquad);
+         for (int iq = 0; iq < nquad; ++iq)
+         {
+            for (int ic = 0; ic < nc; ++ic)
+            {
+               xpt(ic) = xquad(iq,ic);
+            }
+            d(xpt, ypt);
+            double w = pa->GetMetricTerms().Weight(i, iq);
+            for (int ic = 0; ic < nc; ++ic)
+            {
+               yquad(iq,ic) = ypt(ic)*w;
+            }
+         }
+         MultAtB(B, yquad, yel);
+         y.AddElementVector(vdofs, yel.Data());
+      }
+   }
+};
+
+template <typename D>
+class GtDB : public Operator
+{
+   const PartialAssembly *pa;
+   const FiniteElementSpace *fes;
+public:
+   D d;
+   GtDB(const PartialAssembly *pa_) : pa(pa_), fes(pa->GetFES()) { }
+   void Mult(const Vector &x, Vector &y) const
+   {
+      const int nc = d.NComponents();
+      const int dim = fes->GetFE(0)->GetDim();
+      Array<int> vdofs;
+      DenseMatrix xquad, xel, yel;
+      DenseTensor yquad;
+      Vector xpt(nc);
+      DenseMatrix F(nc,dim);
+      for (int i = 0; i < fes->GetNE(); i++)
+      {
+         const DenseMatrix &B = pa->BasisEval(i);
+         const DenseTensor &G = pa->DerivEval(i);
+         const int nquad = B.Height();
+         const int ndof = B.Width();
+         xquad.SetSize(nquad, nc);
+         yquad.SetSize(nquad, nc, dim);
+         xel.SetSize(ndof, nc);
+         yel.SetSize(ndof, nc);
+         fes->GetElementVDofs(i, vdofs);
+         x.GetSubVector(vdofs, xel.Data());
+         mfem::Mult(B, xel, xquad);
+         yquad = 0.0;
+         for (int iq = 0; iq < nquad; ++iq)
+         {
+            for (int ic = 0; ic < nc; ++ic)
+            {
+               xpt(ic) = xquad(iq, ic);
+            }
+            d(xpt, F.Data());
+            const DenseMatrix &Jinv = pa->GetMetricTerms().InverseJacobian(i, iq);
+            double w = pa->GetMetricTerms().Weight(i, iq);
+            for (int ic = 0; ic < nc; ++ic)
+            {
+               for (int d2 = 0; d2 < dim; ++d2)
+               {
+                  for (int d1 = 0; d1 < dim; ++d1)
+                  {
+                     yquad(iq, ic, d1) += w*Jinv(d1,d2)*F(ic,d2);
+                  }
+               }
+            }
+         }
+         for (int d = 0; d < dim; ++d)
+         {
+            MultAtB(G(d), yquad(d), yel);
+            y.AddElementVector(vdofs, yel.Data());
+         }
+      }
+   }
+};
+
+template <typename D>
+class BtDB_face : public Operator
+{
+   const PartialAssembly *pa;
+   const FiniteElementSpace *fes;
+public:
+   D d;
+   BtDB_face(const PartialAssembly *pa_) : pa(pa_), fes(pa->GetFES()) { }
+   void Mult(const Vector &x, Vector &y) const
+   {
+      const int nfaces = fes->GetMesh()->GetNumFaces();
+      const int nc = d.NComponents();
+      Array<int> vdofs1, vdofs2;
+      DenseMatrix xel1, xel2, xquad1, xquad2, yel1, yel2, yquad;
+      Vector xpt1(nc), xpt2(nc), Fdotn(nc);
+      for (int i = 0; i < nfaces; ++i)
+      {
+         PartialAssembly::F2E i1, i2;
+         int iel1, iel2;
+         pa->GetFES()->GetMesh()->GetFaceElements(i, &iel1, &iel2);
+         pa->GetF2Es(i, i1, i2);
+
+         // 1. Evaluate DOFs from element 1 at the face
+         const DenseMatrix &Bface1 = pa->FaceEval(i1);
+         const int nquad = Bface1.Height();
+         const int ndof1 = Bface1.Width();
+         xel1.SetSize(ndof1,nc);
+         yel1.SetSize(ndof1,nc);
+         xquad1.SetSize(nquad,nc);
+         fes->GetElementVDofs(iel1, vdofs1);
+         x.GetSubVector(vdofs1, xel1.Data());
+         mfem::Mult(Bface1, xel1, xquad1);
+         // Is the face interior? This will check validity of the second element
+         if (i2)
+         {
+            // 2. If the face is interior, evaluate DOFs from element 2 at the face
+            const DenseMatrix &Bface2 = pa->FaceEval(i2);
+            const int ndof2 = Bface2.Width();
+            xel2.SetSize(ndof2, nc);
+            yel2.SetSize(ndof2, nc);
+            fes->GetElementVDofs(iel2, vdofs2);
+            x.GetSubVector(vdofs2, xel2.Data());
+            xquad2.SetSize(nquad, nc);
+            mfem::Mult(Bface2, xel2, xquad2);
+         }
+         else
+         {
+            xquad2 = xquad1;
+         }
+
+         // 3. Evaluate the two-point numerical flux at each quadrature point
+         yquad.SetSize(nquad, nc);
+         for (int iq = 0; iq < nquad; ++iq)
+         {
+            for (int ic = 0; ic < nc; ++ic)
+            {
+               xpt1(ic) = xquad1(iq, ic);
+               xpt2(ic) = xquad2(iq, ic);
+            }
+            // Compute face normals
+            const Vector &nvec = pa->GetMetricTerms().Normal(i, iq);
+            d(xpt1, xpt2, nvec.GetData(), Fdotn);
+            double w = pa->GetMetricTerms().FaceWeight(i, iq);
+            // Multiply by geometric factors and quadrature weights
+            for (int ic = 0; ic < nc; ++ic)
+            {
+               yquad(iq, ic) = w*Fdotn(ic);
+            }
+         }
+
+         // 4. Integrate against test functions by multiplying pointwise by
+         //    quadrature weights, and then multiply by the transpose of the
+         //    operators from steps 1 and 2
+         // 5. Add back to residual DOFs for element 1 (and element 2
+         //    if the face is interior)
+
+         MultAtB(Bface1, yquad, yel1);
+         yel1 *= -1.0;
+         y.AddElementVector(vdofs1, yel1.Data());
+
+         if (i2)
+         {
+            const DenseMatrix &Bface2 = pa->FaceEval(i2);
+            MultAtB(Bface2, yquad, yel2);
+            y.AddElementVector(vdofs2, yel2.Data());
+         }
+      }
+   }
+};
+
+} // namespace dg
+} // namespace mfem
+
+#endif

miniapps/conslaw/include/dg_pa.hpp

@@ -0,0 +1,74 @@
+#ifndef DGPA_PA
+#define DGPA_PA
+
+#include "mfem.hpp"
+#include "dg_metric.hpp"
+#include <unordered_map>
+
+namespace mfem
+{
+namespace dg
+{
+
+class PartialAssembly
+{
+public:
+   struct F2E
+   {
+      int face_type;
+      int elem_type;
+      int info;
+      operator bool() const { return face_type != -1; }
+   };
+private:
+   struct F2EHash
+   {
+      std::size_t operator()(const F2E &i) const
+      {
+         static std::hash<int> h;
+         return h(i.face_type)^h(i.elem_type)^h(i.info);
+      }
+   };
+   struct F2EEq
+   {
+      bool operator()(const F2E &i1, const F2E &i2) const
+      {
+         return (i1.face_type == i2.face_type) && (i1.elem_type == i2.elem_type)
+                && (i1.info == i2.info);
+      }
+   };
+
+
+   const FiniteElementSpace *fes;
+   IntegrationRule ir, ir_face;
+   MetricTerms metric;
+   DenseMatrix B;
+   DenseTensor G;
+   std::unordered_map<F2E,DenseMatrix,F2EHash,F2EEq> Bfaces;
+
+   void FormFaceEvaluation(const FiniteElement *fe,
+                           IntegrationPointTransformation *loc,
+                           DenseMatrix &Bface);
+public:
+   PartialAssembly(const FiniteElementSpace *fes_);
+   const FiniteElementSpace *GetFES() const;
+
+   int NQuad(int iel) const;
+
+   // Partial assembly operators
+   const DenseMatrix& BasisEval(int iel) const;
+   const DenseTensor& DerivEval(int iel) const;
+   const DenseMatrix& FaceEval(const F2E &i) const;
+
+   // Face info
+   void GetF2Es(int fid, F2E &i1, F2E &i2) const;
+
+   // Metric terms
+   const MetricTerms& GetMetricTerms() const { return metric; }
+};
+
+
+} // namespace dg
+} // namespace mfem
+
+#endif

miniapps/conslaw/makefile

@@ -0,0 +1,79 @@
+MFEM_DIR ?= ../..
+MFEM_BUILD_DIR ?= ../..
+CONFIG_MK = $(MFEM_BUILD_DIR)/config/config.mk
+
+MFEM_LIB_FILE = mfem_is_not_built
+-include $(CONFIG_MK)
+
+# Directory structure
+LIB_DIR=./lib
+SOURCE_DIR=src
+INC_DIR=include
+BUILD_DIR=build
+APP_DIR=apps
+DEPENDS=$(BUILD_DIR)/.depends
+
+# dgpa shared library
+LIB_NAME=dgpa
+LIB=$(LIB_DIR)/lib$(LIB_NAME).a
+LIB_SOURCES=$(wildcard $(SOURCE_DIR)/*.cpp)
+LIB_OBJECTS=$(LIB_SOURCES:$(SOURCE_DIR)/%.cpp=$(BUILD_DIR)/$(LIB_DIR)/%.o)
+
+# "apps" using the dgpa library
+APP_SOURCES=$(wildcard $(APP_DIR)/*.cpp)
+APP_OBJECTS=$(APP_SOURCES:$(APP_DIR)/%.cpp=$(BUILD_DIR)/$(APP_DIR)/%.o)
+APPS=$(APP_SOURCES:%.cpp=%)
+
+# Compiler configuration
+CXXFLAGS= -std=c++11 ${MFEM_CXXFLAGS}
+LFLGAS=${MFEM_LIBS} ${MFEM_EXT_LIBS}
+INCFLAGS=-I$(SOURCE_DIR) -I${INC_DIR} ${MFEM_INCFLAGS}
+
+.PHONY: all clean style
+all: $(LIB) $(APPS)
+
+# Build all the apps (the library is a dependency)
+$(APPS): $(APP_DIR)/%: $(BUILD_DIR)/$(APP_DIR)/%.o $(LIB)
+	$(MFEM_CXX) $< $(LFLGAS) -L$(LIB_DIR) -l$(LIB_NAME) -o $@
+
+$(BUILD_DIR)/$(APP_DIR)/%.o: makefile | $(BUILD_DIR)
+	$(MFEM_CXX) -c $(CXXFLAGS) $(INCFLAGS) -o $@ $(APP_DIR)/$*.cpp
+
+# Build the library
+$(LIB): $(LIB_OBJECTS)
+	$(AR) $(ARFLAGS) $(@) $(LIB_OBJECTS)
+
+$(BUILD_DIR)/$(LIB_DIR)/%.o: makefile | $(BUILD_DIR)
+	$(MFEM_CXX) -c $(CXXFLAGS) $(INCFLAGS) -o $@ $(SOURCE_DIR)/$*.cpp
+
+# Use the compiler to determine dependencies on header files
+# Some awk magic in the next target
+# Prefix all lines matching the regexp /^.*\.o/ with build/ (to match e.g. "file.o:")
+# and leave other lines alone
+$(DEPENDS): $(APP_SOURCES) $(LIB_SOURCES) | $(BUILD_DIR)
+	$(CXX) -std=c++11 $(INCFLAGS) -MM $(LIB_SOURCES) \
+		| awk '/^.*\.o:/{ print "$(BUILD_DIR)/$(LIB_DIR)/" $$0; next } 1'  >$@
+	$(CXX) -std=c++11 $(INCFLAGS) -MM $(APP_SOURCES) \
+		| awk '/^.*\.o:/{ print "$(BUILD_DIR)/$(APP_DIR)/" $$0; next } 1'  >>$@
+
+# Rebuild dependencies unless doing "make clean" or "make style"
+ifneq ($(MAKECMDGOALS),clean)
+ifneq ($(MAKECMDGOALS),style)
+-include $(DEPENDS)
+endif
+endif
+
+$(BUILD_DIR):
+	mkdir -p $@/$(APP_DIR)
+	mkdir -p $@/$(LIB_DIR)
+	mkdir -p $(LIB_DIR)
+
+clean:
+	rm -rf $(APPS) $(LIB) $(BUILD_DIR)
+
+FORMAT_FILES = $(foreach dir,$(SOURCE_DIR) $(INC_DIR) $(APP_DIR),"$(dir)/*.?pp")
+ASTYLE = astyle --options=$(MFEM_DIR)/config/mfem.astylerc
+style:
+	@if ! $(ASTYLE) $(FORMAT_FILES) | grep Formatted; then\
+	   echo "No source files were changed.";\
+	fi

miniapps/conslaw/src/dg_mass.cpp

@@ -0,0 +1,88 @@
+#include "dg_mass.hpp"
+
+namespace mfem
+{
+namespace dg
+{
+
+Mass::Mass(const PartialAssembly *pa_)
+   : pa(pa_),
+     fes(pa->GetFES())
+{
+   // Precompute mass matrix
+   // For now assume constant number of DOFs per element (i.e. uniform p)
+   const int ndof = fes->GetFE(0)->GetDof();
+   M.SetSize(ndof, ndof, fes->GetNE());
+
+   // Assemble each local mass matrix and insert into the dense tensor
+   MassIntegrator mi;
+   for (int i = 0; i < fes->GetNE(); i++)
+   {
+      mi.AssembleElementMatrix(*fes->GetFE(i),
+                               *fes->GetElementTransformation(i),
+                               M(i));
+   }
+}
+
+const PartialAssembly* Mass::GetPA() const
+{
+   return pa;
+}
+
+void Mass::Mult(const Vector &x, Vector &y) const
+{
+   Array<int> vdofs;
+   const int ndof = fes->GetFE(0)->GetDof();
+   DenseMatrix xel, yel;
+   for (int i = 0; i < fes->GetNE(); i++)
+   {
+      fes->GetElementVDofs(i, vdofs);
+      int nc = vdofs.Size()/ndof;
+      xel.SetSize(ndof, nc);
+      yel.SetSize(ndof, nc);
+      x.GetSubVector(vdofs, xel.Data());
+      mfem::Mult(M(i), xel, yel);
+      y.SetSubVector(vdofs, yel.Data());
+   }
+}
+
+MassInverse::MassInverse(const Mass *mass_)
+   : mass(mass_),
+     fes(mass->fes)
+{
+   // For now assume constant number of DOFs per element (i.e. uniform p)
+   const int ndof = fes->GetFE(0)->GetDof();
+   Minv.SetSize(ndof, ndof, fes->GetNE());
+
+   // Extract the local mass matrices and then invert, inserting the
+   // result into the Minv dense tensor
+   DenseMatrix Me(ndof);
+   DenseMatrixInverse Me_inv(&Me);
+   for (int i = 0; i < fes->GetNE(); i++)
+   {
+      Me = mass->M(i);
+      Me_inv.Factor();
+      Me_inv.GetInverseMatrix(Minv(i));
+   }
+}
+
+void MassInverse::Mult(const Vector &x, Vector &y) const
+{
+   // y must already be allocated/sized here (is that correct?)
+   Array<int> vdofs;
+   const int ndof = fes->GetFE(0)->GetDof();
+   DenseMatrix xel, yel;
+   for (int i = 0; i < fes->GetNE(); i++)
+   {
+      fes->GetElementVDofs(i, vdofs);
+      int nc = vdofs.Size()/ndof;
+      xel.SetSize(ndof, nc);
+      yel.SetSize(ndof, nc);
+      x.GetSubVector(vdofs, xel.Data());
+      mfem::Mult(Minv(i), xel, yel);
+      y.SetSubVector(vdofs, yel.Data());
+   }
+}
+
+} // namespace dg
+} // namespace mfem

miniapps/conslaw/src/dg_metric.cpp

@@ -0,0 +1,109 @@
+#include "dg_metric.hpp"
+
+namespace mfem
+{
+namespace dg
+{
+
+MetricTerms::MetricTerms()
+{
+   nel = 0;
+   nquad = 0;
+   nfaces = 0;
+   dim = 0;
+   nquad_face = 0;
+}
+
+void MetricTerms::Precompute(const FiniteElementSpace *fes,
+                             const IntegrationRule *ir,
+                             const IntegrationRule *ir_face)
+{
+   Mesh *mesh = fes->GetMesh();
+
+   nquad = ir->Size();
+   nquad_face = ir_face->Size();
+   nel = fes->GetNE();
+   dim = fes->GetFE(0)->GetDim();
+   nfaces = mesh->GetNumFaces();
+
+   w.SetSize(nquad);
+   for (int iq = 0; iq < nquad; ++iq)
+   {
+      w(iq) = ir->IntPoint(iq).weight;
+   }
+
+   detJ.SetSize(nel*nquad);
+   Jinv.SetSize(dim, dim, nel*nquad);
+
+   for (int i = 0; i < nel; ++i)
+   {
+      ElementTransformation *tr = fes->GetElementTransformation(i);
+      for (int iq = 0; iq < nquad; ++iq)
+      {
+         int idx = i*nquad + iq;
+         const IntegrationPoint &ip = ir->IntPoint(iq);
+         tr->SetIntPoint(&ip);
+         detJ(idx) = tr->Weight();
+         Jinv(idx) = tr->InverseJacobian();
+      }
+   }
+
+   wface.SetSize(nquad_face);
+   for (int iq = 0; iq < nquad_face; ++iq)
+   {
+      wface(iq) = ir_face->IntPoint(iq).weight;
+   }
+
+   nvec.SetSize(dim, nfaces*nquad_face);
+   for (int i = 0; i < nfaces; ++i)
+   {
+      FaceElementTransformations *tr = mesh->GetFaceElementTransformations(i);
+      for (int iq = 0; iq < nquad_face; ++iq)
+      {
+         int idx = i*nquad_face + iq;
+         const IntegrationPoint &ip = ir_face->IntPoint(iq);
+         tr->Face->SetIntPoint(&ip);
+         if (dim == 1)
+         {
+            IntegrationPoint eip1;
+            tr->Loc1.Transform(ip, eip1);
+            nvec(0,idx) = 2*eip1.x - 1.0;
+         }
+         else
+         {
+            Vector normal;
+            nvec.GetColumnReference(idx, normal);
+            CalcOrtho(tr->Face->Jacobian(), normal);
+         }
+      }
+   }
+}
+
+double MetricTerms::JacobianDeterminant(int elid, int iq) const
+{
+   return detJ(elid*nquad + iq);
+}
+
+double MetricTerms::Weight(int elid, int iq) const
+{
+   return w(iq)*JacobianDeterminant(elid, iq);
+}
+
+const DenseMatrix& MetricTerms::InverseJacobian(int elid, int iq) const
+{
+   return Jinv(elid*nquad + iq);
+}
+
+const Vector& MetricTerms::Normal(int fid, int iq) const
+{
+   nvec.GetColumnReference(fid*nquad_face + iq, n);
+   return n;
+}
+
+double MetricTerms::FaceWeight(int fid, int iq) const
+{
+   return wface(iq);
+}
+
+} // namespace dg
+} // namespace mfem

miniapps/conslaw/src/dg_pa.cpp

@@ -0,0 +1,136 @@
+#include "dg_pa.hpp"
+
+namespace mfem
+{
+namespace dg
+{
+
+void PartialAssembly::FormFaceEvaluation(const FiniteElement *fe,
+                                         IntegrationPointTransformation *loc,
+                                         DenseMatrix &Bface)
+{
+   const int nquad = ir_face.Size();
+   const int ndof = fe->GetDof();
+   Vector shape(ndof);
+   Bface.SetSize(nquad, ndof);
+   for (int iq = 0; iq < nquad; ++iq)
+   {
+      IntegrationPoint ip;
+      loc->Transform(ir_face.IntPoint(iq), ip);
+      fe->CalcShape(ip, shape);
+      for (int i = 0; i < ndof; ++i)
+      {
+         Bface(iq,i) = shape(i);
+      }
+   }
+}
+
+PartialAssembly::PartialAssembly(const FiniteElementSpace *fes_)
+   : fes(fes_)
+{
+   int order = 3*fes->GetFE(0)->GetOrder();
+   ir = IntRules.Get(fes->GetFE(0)->GetGeomType(), order);
+   ir_face = IntRules.Get(fes->GetMesh()->GetFaceGeometryType(0), order);
+   metric.Precompute(fes, &ir, &ir_face);
+
+   // Precompute basis evaluation matrix B and basis derivative matrix G
+   const int dim = fes->GetFE(0)->GetDim();
+   const int nquad = ir.Size();
+   const int ndof = fes->GetFE(0)->GetDof();
+   B.SetSize(nquad, ndof);
+   G.SetSize(nquad, ndof, dim);
+
+   const FiniteElement *fe = fes->GetFE(0);
+   Vector shape(ndof);
+   DenseMatrix dshape(ndof, dim);
+   for (int iq = 0; iq < nquad; ++iq)
+   {
+      const IntegrationPoint &ip = ir.IntPoint(iq);
+      fe->CalcShape(ip, shape);
+      fe->CalcDShape(ip, dshape);
+      for (int i = 0; i < ndof; ++i)
+      {
+         B(iq,i) = shape(i);
+         for (int d = 0; d < dim; ++ d)
+         {
+            G(iq, i, d) = dshape(i, d);
+         }
+      }
+   }
+
+   // Precompute face evaluation matrices
+   Mesh *mesh = fes->GetMesh();
+   const int nfaces = mesh->GetNumFaces();
+   for (int i = 0; i < nfaces; ++i)
+   {
+      DenseMatrix Bface1, Bface2;
+      F2E i1, i2;
+
+      FaceElementTransformations *tr = mesh->GetFaceElementTransformations(i);
+      GetF2Es(i, i1, i2);
+
+      if (Bfaces.find(i1) == Bfaces.end())
+      {
+         FormFaceEvaluation(fes->GetFE(tr->Elem1No), &tr->Loc1, Bface1);
+         Bfaces[i1] = Bface1;
+      }
+
+      if (tr->Elem2No >= 0)
+      {
+         if (Bfaces.find(i2) == Bfaces.end())
+         {
+            FormFaceEvaluation(fes->GetFE(tr->Elem2No), &tr->Loc2, Bface2);
+            Bfaces[i2] = Bface2;
+         }
+      }
+   }
+}
+
+int PartialAssembly::NQuad(int iel) const
+{
+   return ir.Size();
+}
+
+const FiniteElementSpace* PartialAssembly::GetFES() const
+{
+   return fes;
+}
+
+void PartialAssembly::GetF2Es(int fid, F2E &i1, F2E &i2) const
+{
+   Mesh *mesh = fes->GetMesh();
+   int iel1, iel2;
+   mesh->GetFaceElements(fid, &iel1, &iel2);
+   i1.elem_type = mesh->GetElementType(iel1);
+   i1.face_type = mesh->GetFaceElementType(fid);
+   mesh->GetFaceInfos(fid, &i1.info, &i2.info);
+   if (iel2 >= 0)
+   {
+      i2.elem_type = mesh->GetElementType(iel2);
+      i2.face_type = i1.face_type;
+   }
+   else
+   {
+      // Invalid
+      i2.info = -1;
+      i2.elem_type = -1;
+      i2.face_type = -1;
+   }
+}
+
+const DenseMatrix& PartialAssembly::BasisEval(int iel) const
+{
+   return B;
+}
+const DenseTensor& PartialAssembly::DerivEval(int iel) const
+{
+   return G;
+}
+
+const DenseMatrix& PartialAssembly::FaceEval(const F2E &i) const
+{
+   return Bfaces.at(i);
+}
+
+} // namespace dg
+} // namespace mfem