Added support for variable-order spaces.

data/circle_3_4.mesh

@@ -0,0 +1,324 @@
+MFEM mesh v1.0
+
+#
+# MFEM Geometry Types (see mesh/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+# PRISM       = 6
+#
+
+dimension
+2
+
+elements
+15
+1 3 3 6 17 9
+1 3 17 7 4 8
+1 3 9 17 8 5
+1 3 0 10 18 11
+1 3 11 18 6 3
+1 3 12 1 13 19
+1 3 19 13 4 7
+1 3 2 14 20 15
+1 3 14 5 8 20
+1 3 20 8 4 13
+1 3 15 20 13 1
+1 3 0 11 21 16
+1 3 11 3 9 21
+1 3 21 9 5 14
+1 3 16 21 14 2
+
+boundary
+12
+1 1 6 17
+2 1 17 7
+3 1 0 10
+4 1 10 18
+5 1 18 6
+6 1 19 12
+7 1 12 1
+8 1 7 19
+9 1 15 2
+10 1 1 15
+11 1 16 0
+12 1 2 16
+
+vertices
+22
+
+nodes
+FiniteElementSpace
+FiniteElementCollection: H1_2D_P4
+VDim: 2
+Ordering: 1
+
+-0.707106781186545 -0.707106781186545
+0.707106781186545 0.707106781186545
+-0.707106781186545 0.707106781186545
+-0.3535533905932725 -0.3535533905932725
+0.3535533905932725 0.3535533905932725
+-0.3535533905932725 0.3535533905932725
+0 -0.3535533905932725
+0.3535533905932725 0
+0 0.3535533905932725
+-0.3535533905932725 0
+0 -1
+-0.5303300858899087 -0.5303300858899087
+1 0
+0.5303300858899087 0.5303300858899087
+-0.5303300858899087 0.5303300858899087
+0 1
+-1 0
+0 0
+0 -0.6721895948517682
+0.6721895948517682 0
+-2.45326946669338e-17 0.6721895948517682
+-0.6721895948517682 -2.45326946669338e-17
+-0.2925042077682053 -0.3535533905932725
+-0.1767766952966363 -0.3535533905932725
+-0.06104918282506723 -0.3535533905932725
+0 -0.2925042077682053
+0 -0.1767766952966363
+0 -0.06104918282506723
+-0.2925042077682053 -1.937718899096505e-18
+-0.1767766952966363 -2.168768018011117e-18
+-0.06104918282506723 -6.942024631559806e-19
+-0.3535533905932725 -0.2925042077682053
+-0.3535533905932725 -0.1767766952966363
+-0.3535533905932725 -0.06104918282506723
+0.2925042077682053 -2.25773586854044e-19
+0.1767766952966363 2.904835825740127e-19
+0.06104918282506723 4.044254824131427e-19
+0.3535533905932725 0.2925042077682053
+0.3535533905932725 0.1767766952966363
+0.3535533905932725 0.06104918282506723
+0.2925042077682053 0.3535533905932725
+0.1767766952966363 0.3535533905932725
+0.06104918282506723 0.3535533905932725
+0 0.2925042077682053
+0 0.1767766952966363
+0 0.06104918282506723
+-0.2925042077682053 0.3535533905932725
+-0.1767766952966363 0.3535533905932725
+-0.06104918282506723 0.3535533905932725
+-0.3535533905932725 0.2925042077682053
+-0.3535533905932725 0.1767766952966363
+-0.3535533905932725 0.06104918282506723
+-0.6145094428537778 -0.789413761659266
+-0.3959323428473652 -0.9179997269052453
+-0.1410304266231549 -0.989913000881131
+0 -0.9426395767133444
+0 -0.8348063668130653
+0 -0.7280766470976995
+-0.4493360908158245 -0.5727497127416228
+-0.2803631562583237 -0.6346920242481003
+-0.09836326743220743 -0.6676463098413217
+-0.6765821897740132 -0.6765821897740132
+-0.6187184335382269 -0.6187184335382269
+-0.5608546773024407 -0.5608546773024407
+0 -0.4080066739179119
+0 -0.51185560644339
+0 -0.6165744197248449
+-0.3840779820058043 -0.3840779820058043
+-0.4419417382415907 -0.4419417382415907
+-0.4998054944773769 -0.4998054944773769
+0.789413761659266 0.6145094428537778
+0.9179997269052453 0.3959323428473652
+0.989913000881131 0.1410304266231549
+0.6765821897740132 0.6765821897740132
+0.6187184335382269 0.6187184335382269
+0.5608546773024407 0.5608546773024407
+0.5727497127416228 0.4493360908158245
+0.6346920242481003 0.2803631562583237
+0.6676463098413217 0.09836326743220743
+0.9426395767133444 0
+0.8348063668130653 0
+0.7280766470976995 0
+0.3840779820058043 0.3840779820058043
+0.4419417382415907 0.4419417382415907
+0.4998054944773769 0.4998054944773769
+0.4080066739179119 0
+0.51185560644339 0
+0.6165744197248449 0
+-0.6765821897740132 0.6765821897740132
+-0.6187184335382269 0.6187184335382269
+-0.5608546773024407 0.5608546773024407
+-0.4493360908158245 0.5727497127416228
+-0.2803631562583237 0.6346920242481003
+-0.09836326743220743 0.6676463098413217
+1.38261702840016e-17 0.9426395767133444
+2.635181002135987e-18 0.8348063668130653
+-1.884789714682629e-17 0.7280766470976995
+-0.6145094428537778 0.789413761659266
+-0.3959323428473652 0.9179997269052453
+-0.1410304266231549 0.989913000881131
+-0.3840779820058043 0.3840779820058043
+-0.4419417382415907 0.4419417382415907
+-0.4998054944773769 0.4998054944773769
+5.694430050849992e-18 0.4080066739179119
+-9.046264100643195e-18 0.51185560644339
+-2.406637988827979e-17 0.6165744197248449
+0.4493360908158245 0.5727497127416228
+0.2803631562583237 0.6346920242481003
+0.09836326743220743 0.6676463098413217
+0.6145094428537778 0.789413761659266
+0.3959323428473652 0.9179997269052453
+0.1410304266231549 0.989913000881131
+-0.5727497127416228 -0.4493360908158245
+-0.6346920242481003 -0.2803631562583237
+-0.6676463098413217 -0.09836326743220743
+-0.9426395767133444 1.38261702840016e-17
+-0.8348063668130653 2.635181002135987e-18
+-0.7280766470976995 -1.884789714682629e-17
+-0.789413761659266 -0.6145094428537778
+-0.9179997269052453 -0.3959323428473652
+-0.989913000881131 -0.1410304266231549
+-0.4080066739179119 5.694430050849992e-18
+-0.51185560644339 -9.046264100643195e-18
+-0.6165744197248449 -2.406637988827979e-17
+-0.5727497127416228 0.4493360908158245
+-0.6346920242481003 0.2803631562583237
+-0.6676463098413217 0.09836326743220743
+-0.789413761659266 0.6145094428537778
+-0.9179997269052453 0.3959323428473652
+-0.989913000881131 0.1410304266231549
+-0.2925042077682053 -0.2925042077682053
+-0.1767766952966363 -0.2925042077682053
+-0.06104918282506723 -0.2925042077682053
+-0.2925042077682053 -0.1767766952966363
+-0.1767766952966363 -0.1767766952966363
+-0.06104918282506723 -0.1767766952966363
+-0.2925042077682053 -0.06104918282506723
+-0.1767766952966363 -0.06104918282506723
+-0.06104918282506723 -0.06104918282506723
+0.06104918282506723 0.06104918282506723
+0.1767766952966363 0.06104918282506723
+0.2925042077682053 0.06104918282506723
+0.06104918282506723 0.1767766952966363
+0.1767766952966363 0.1767766952966363
+0.2925042077682053 0.1767766952966363
+0.06104918282506723 0.2925042077682053
+0.1767766952966363 0.2925042077682053
+0.2925042077682053 0.2925042077682053
+-0.2925042077682053 0.06104918282506723
+-0.1767766952966363 0.06104918282506723
+-0.06104918282506723 0.06104918282506723
+-0.2925042077682053 0.1767766952966363
+-0.1767766952966363 0.1767766952966363
+-0.06104918282506723 0.1767766952966363
+-0.2925042077682053 0.2925042077682053
+-0.1767766952966363 0.2925042077682053
+-0.06104918282506723 0.2925042077682053
+-0.5853433152017873 -0.7522063629666647
+-0.3750499116779216 -0.8688976983986201
+-0.1332489509046759 -0.9335917585418947
+-0.530810235275664 -0.681431690503692
+-0.3365495705701402 -0.7760363182652825
+-0.118982888320896 -0.8276309647697756
+-0.4772305101537735 -0.6103572471372083
+-0.2994187471166691 -0.6834400186968074
+-0.1053286339663329 -0.7226539726955523
+-0.4216899203778937 -0.5350670333567459
+-0.2616641810583542 -0.5860082222568944
+-0.09155721275052374 -0.6128808628755988
+-0.3699443285846694 -0.4634369827877818
+-0.2271681353001397 -0.4938868826316494
+-0.07907999957616745 -0.5096954174503119
+-0.3190349108655897 -0.3915624660435183
+-0.1938732741916287 -0.4019674242664178
+-0.06713939478447406 -0.4072844098734093
+0.9335917585418947 0.1332489509046759
+0.8688976983986201 0.3750499116779216
+0.7522063629666647 0.5853433152017873
+0.8276309647697756 0.118982888320896
+0.7760363182652825 0.3365495705701402
+0.681431690503692 0.530810235275664
+0.7226539726955523 0.1053286339663329
+0.6834400186968074 0.2994187471166691
+0.6103572471372083 0.4772305101537735
+0.6128808628755988 0.09155721275052374
+0.5860082222568944 0.2616641810583542
+0.5350670333567459 0.4216899203778937
+0.5096954174503119 0.07907999957616745
+0.4938868826316494 0.2271681353001397
+0.4634369827877818 0.3699443285846694
+0.4072844098734093 0.06713939478447406
+0.4019674242664178 0.1938732741916287
+0.3915624660435183 0.3190349108655861
+-0.5853433152017873 0.7522063629666647
+-0.530810235275664 0.681431690503692
+-0.4772305101537735 0.6103572471372083
+-0.3750499116779216 0.8688976983986201
+-0.3365495705701402 0.7760363182652825
+-0.2994187471166691 0.6834400186968074
+-0.1332489509046759 0.9335917585418947
+-0.118982888320896 0.8276309647697756
+-0.1053286339663329 0.7226539726955523
+-0.4216899203778937 0.5350670333567459
+-0.3699443285846694 0.4634369827877818
+-0.3190349108655861 0.3915624660435183
+-0.2616641810583542 0.5860082222568944
+-0.2271681353001397 0.4938868826316494
+-0.1938732741916287 0.4019674242664178
+-0.09155721275052374 0.6128808628755988
+-0.07907999957616745 0.5096954174503119
+-0.06713939478447406 0.4072844098734093
+0.09155721275052374 0.6128808628755988
+0.07907999957616745 0.5096954174503119
+0.06713939478447406 0.4072844098734093
+0.2616641810583542 0.5860082222568944
+0.2271681353001397 0.4938868826316494
+0.1938732741916287 0.4019674242664178
+0.4216899203778937 0.5350670333567459
+0.3699443285846694 0.4634369827877818
+0.3190349108655897 0.3915624660435183
+0.1332489509046759 0.9335917585418947
+0.118982888320896 0.8276309647697756
+0.1053286339663329 0.7226539726955523
+0.3750499116779216 0.8688976983986201
+0.3365495705701402 0.7760363182652825
+0.2994187471166691 0.6834400186968074
+0.5853433152017873 0.7522063629666647
+0.530810235275664 0.681431690503692
+0.4772305101537735 0.6103572471372083
+-0.7522063629666647 -0.5853433152017873
+-0.681431690503692 -0.530810235275664
+-0.6103572471372083 -0.4772305101537735
+-0.8688976983986201 -0.3750499116779216
+-0.7760363182652825 -0.3365495705701402
+-0.6834400186968074 -0.2994187471166691
+-0.9335917585418947 -0.1332489509046759
+-0.8276309647697756 -0.118982888320896
+-0.7226539726955523 -0.1053286339663329
+-0.5350670333567459 -0.4216899203778937
+-0.4634369827877818 -0.3699443285846694
+-0.3915624660435183 -0.3190349108655861
+-0.5860082222568944 -0.2616641810583542
+-0.4938868826316494 -0.2271681353001397
+-0.4019674242664178 -0.1938732741916287
+-0.6128808628755988 -0.09155721275052374
+-0.5096954174503119 -0.07907999957616745
+-0.4072844098734093 -0.06713939478447406
+-0.6128808628755988 0.09155721275052374
+-0.5096954174503119 0.07907999957616745
+-0.4072844098734093 0.06713939478447406
+-0.5860082222568944 0.2616641810583542
+-0.4938868826316494 0.2271681353001397
+-0.4019674242664178 0.1938732741916287
+-0.5350670333567459 0.4216899203778937
+-0.4634369827877818 0.3699443285846694
+-0.3915624660435183 0.3190349108655897
+-0.9335917585418947 0.1332489509046759
+-0.8276309647697756 0.118982888320896
+-0.7226539726955523 0.1053286339663329
+-0.8688976983986201 0.3750499116779216
+-0.7760363182652825 0.3365495705701402
+-0.6834400186968074 0.2994187471166691
+-0.7522063629666647 0.5853433152017873
+-0.681431690503692 0.530810235275664
+-0.6103572471372083 0.4772305101537735

examples/advection.cpp

@@ -0,0 +1,367 @@
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+#include <algorithm>
+
+using namespace std;
+using namespace mfem;
+
+
+void Prefine(FiniteElementSpace & fes_old,
+             GridFunction &u, Coefficient &gf_ex, GridFunction &orders_gf,
+             double min_thresh, double max_thresh);
+
+
+// Velocity coefficient
+void velocity_function(const Vector &x, Vector &v);
+
+// Initial condition
+double u0_function(const Vector &x, double);
+
+// Inflow boundary condition
+double inflow_function(const Vector &x);
+
+// Mesh bounding box
+Vector bb_min, bb_max;
+
+class FE_Evolution : public TimeDependentOperator
+{
+private:
+   BilinearForm &M, &K;
+   const Vector &b;
+   Solver *M_prec;
+   CGSolver M_solver;
+
+   mutable Vector z;
+
+public:
+   FE_Evolution(BilinearForm &M_, BilinearForm &K_, const Vector &b_);
+   void Update();
+   virtual void Mult(const Vector &x, Vector &y) const;
+
+   virtual ~FE_Evolution();
+};
+
+
+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 = 1;
+   double t_final = 10.0;
+   double dt = 0.0005;
+   bool visualization = true;
+   int vis_steps = 5;
+
+   int precision = 8;
+   cout.precision(precision);
+
+   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.");
+   args.AddOption(&order, "-o", "--order",
+                  "Order (degree) of the finite elements.");
+   args.AddOption(&t_final, "-tf", "--t-final",
+                  "Final time; start time is 0.");
+   args.AddOption(&dt, "-dt", "--time-step",
+                  "Time step.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&vis_steps, "-vs", "--visualization-steps",
+                  "Visualize every n-th timestep.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   Mesh mesh0 = Mesh::MakeCartesian2D(64, 1, mfem::Element::QUADRILATERAL,false, 2,
+                                      1);
+
+
+   std::vector<Vector> translations = {Vector({2.0,0.0}), };
+
+
+   Mesh mesh = Mesh::MakePeriodic(mesh0,
+                                  mesh0.CreatePeriodicVertexMapping(translations));
+
+
+   mesh.EnsureNCMesh();
+
+   int dim = mesh.Dimension();
+
+   ODESolver *ode_solver = new RK4Solver;
+
+
+   mesh.GetBoundingBox(bb_min, bb_max, max(order, 1));
+
+   // 5. Define the discontinuous DG finite element space of the given
+   //    polynomial order on the refined mesh.
+   DG_FECollection fec(order, dim, BasisType::GaussLobatto);
+   FiniteElementSpace fes(&mesh, &fec);
+   FiniteElementSpace fes_old(&mesh, &fec);
+
+   cout << "Number of unknowns: " << fes.GetVSize() << endl;
+
+   VectorFunctionCoefficient velocity(dim, velocity_function);
+   FunctionCoefficient inflow(inflow_function);
+   FunctionCoefficient u0(u0_function);
+
+   BilinearForm m(&fes);
+   BilinearForm k(&fes);
+   m.AddDomainIntegrator(new MassIntegrator);
+   constexpr double alpha = -1.0;
+   k.AddDomainIntegrator(new ConvectionIntegrator(velocity, alpha));
+   k.AddInteriorFaceIntegrator(
+      new NonconservativeDGTraceIntegrator(velocity, alpha));
+   k.AddBdrFaceIntegrator(
+      new NonconservativeDGTraceIntegrator(velocity, alpha));
+
+   LinearForm b(&fes);
+   b.AddBdrFaceIntegrator(
+      new BoundaryFlowIntegrator(inflow, velocity, alpha,-0.5));
+
+   m.Assemble();
+   int skip_zeros = 0;
+   k.Assemble(skip_zeros);
+   b.Assemble();
+   m.Finalize();
+   k.Finalize(skip_zeros);
+
+   // 7. Define the initial conditions, save the corresponding grid function to
+   //    a file and (optionally) save data in the VisIt format and initialize
+   //    GLVis visualization.
+   GridFunction u(&fes);
+   u0.SetTime(0.);
+   u.ProjectCoefficient(u0);
+
+   L2_FECollection orders_fec(0,dim);
+   FiniteElementSpace orders_fes(&mesh,&orders_fec);
+   GridFunction orders_gf(&orders_fes);
+   for (int i = 0; i<mesh.GetNE(); i++) { orders_gf(i) = order; }
+
+   socketstream sout;
+   socketstream meshout;
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      sout.open(vishost, visport);
+      meshout.open(vishost, visport);
+      if (!sout)
+      {
+         cout << "Unable to connect to GLVis server at "
+              << vishost << ':' << visport << endl;
+         visualization = false;
+         cout << "GLVis visualization disabled.\n";
+      }
+      else
+      {
+         sout.precision(precision);
+         sout << "solution\n" << mesh << u;
+         sout << flush;
+         meshout.precision(precision);
+         meshout << "solution\n" << mesh << orders_gf;
+         meshout << flush;
+      }
+   }
+
+
+   // 8. Define the time-dependent evolution operator describing the ODE
+   //    right-hand side, and perform time-integration (looping over the time
+   //    iterations, ti, with a time-step dt).
+   FE_Evolution adv(m, k, b);
+
+   double t = 0.0;
+   adv.SetTime(t);
+
+   GridFunction gf_ex(&fes);
+   FunctionCoefficient u_ex(u0_function);
+
+   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;
+         u_ex.SetTime(t);
+         Prefine(fes_old,u,u_ex, orders_gf, 5e-5, 5e-4);
+         m.Update();
+         m.Assemble();
+         m.Finalize();
+         k.Update();
+         k.Assemble(skip_zeros);
+         k.Finalize(skip_zeros);
+         b.Update();
+         b.Assemble();
+         adv.Update();
+         ode_solver->Init(adv);
+         if (visualization)
+         {
+            GridFunction * pr_u = ProlongToMaxOrder(&u);
+            sout << "solution\n" << mesh << *pr_u << flush;
+            meshout << "solution\n" << mesh << orders_gf << flush;
+         }
+      }
+   }
+
+   // 10. Free the used memory.
+   delete ode_solver;
+
+   return 0;
+}
+
+
+// Implementation of class FE_Evolution
+FE_Evolution::FE_Evolution(BilinearForm &M_, BilinearForm &K_, const Vector &b_)
+   : TimeDependentOperator(M_.Height()), M(M_), K(K_), b(b_), z(M_.Height())
+{
+   Array<int> ess_tdof_list;
+   M_prec = new OperatorJacobiSmoother(M, ess_tdof_list);
+   M_solver.SetOperator(M);
+   M_solver.SetPreconditioner(*M_prec);
+   M_solver.iterative_mode = false;
+   M_solver.SetRelTol(1e-9);
+   M_solver.SetAbsTol(0.0);
+   M_solver.SetMaxIter(100);
+   M_solver.SetPrintLevel(0);
+}
+
+void FE_Evolution::Update()
+{
+   height = M.Height();
+   width = M.Width();
+   z.SetSize(M.Height());
+
+   Array<int> ess_tdof_list;
+   delete M_prec;
+   M_prec = new OperatorJacobiSmoother(M, ess_tdof_list);
+   M_solver.SetOperator(M);
+   M_solver.SetPreconditioner(*M_prec);
+   M_solver.iterative_mode = false;
+   M_solver.SetRelTol(1e-9);
+   M_solver.SetAbsTol(0.0);
+   M_solver.SetMaxIter(100);
+   M_solver.SetPrintLevel(0);
+}
+
+void FE_Evolution::Mult(const Vector &x, Vector &y) const
+{
+   // y = M^{-1} (K x + b)
+   K.Mult(x, z);
+   z += b;
+   M_solver.Mult(z, y);
+}
+
+FE_Evolution::~FE_Evolution()
+{
+   delete M_prec;
+}
+
+// Velocity coefficient
+void velocity_function(const Vector &x, Vector &v)
+{
+   v.SetSize(2);
+   v(0) = 1.;
+   v(1) = 0.;
+}
+
+// Initial condition
+double u0_function(const Vector &x, double t)
+{
+   // give x0, y0;
+   double x0 = 0.5;
+   // double y0 = 0.5;
+   double w = 100.;
+   double c = 1.;
+   double ds = c*t;
+
+   double xx = x(0) - ds;
+   double yy = x(1) - ds;
+
+   double tol = 1e-6;
+   if (xx>= 2.0+tol || xx<= 0.0-tol)
+   {
+      xx -= (int)xx;
+   }
+   if (yy>= 1.0+tol || yy<= 0.0-tol)
+   {
+      yy -= (int)yy;
+   }
+
+   double dr2 = (xx-x0)*(xx-x0);
+   return 1. + exp(-w*dr2);
+}
+
+// Inflow boundary condition (zero for the problems considered in this example)
+double inflow_function(const Vector &x)
+{
+   return 1.0;
+}
+
+
+void Prefine(FiniteElementSpace & fes_old,
+             GridFunction &u, Coefficient &ex, GridFunction &orders_gf,
+             double min_thresh, double max_thresh)
+{
+   // get element errors
+   FiniteElementSpace * fes = u.FESpace();
+   int ne = fes->GetMesh()->GetNE();
+   Vector errors(ne);
+   u.ComputeElementL2Errors(ex,errors);
+   for (int i = 0; i<ne; i++)
+   {
+      double error = errors(i);
+      int order = fes->GetElementOrder(i);
+      if (error < min_thresh && order > 1)
+      {
+         fes->SetElementOrder(i,order-1);
+      }
+      else if (error > max_thresh && order < 2)
+      {
+         fes->SetElementOrder(i, order+1);
+      }
+      else
+      {
+         // do nothing
+      }
+   }
+
+   fes->Update(false);
+
+   PRefinementTransferOperator * T = new PRefinementTransferOperator(fes_old,*fes);
+
+   GridFunction u_fine(fes);
+   T->Mult(u,u_fine);
+
+   // copy the orders to the old space
+   for (int i = 0; i<ne; i++)
+   {
+      int order = fes->GetElementOrder(i);
+      fes_old.SetElementOrder(i,order);
+      orders_gf(i) = order;
+   }
+   fes_old.Update(false);
+
+   delete T;
+
+   // update old gridfuntion;
+   u = u_fine;
+
+}

examples/advection_href.cpp

@@ -0,0 +1,706 @@
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+#include <algorithm>
+
+using namespace std;
+using namespace mfem;
+
+double sx, sy;
+
+Vector vel;
+
+void Prefine(FiniteElementSpace & fes_old,
+             GridFunction &u, Coefficient &gf_ex, GridFunction &orders_gf,
+             double min_thresh, double max_thresh);
+
+void Hrefine(GridFunction &u, Coefficient &gf_ex, double min_thresh,
+             double max_thresh);
+
+void Hrefine2(GridFunction &u, Coefficient &gf_ex, double min_thresh,
+              double max_thresh);
+
+Table * Refine(Array<int> ref_actions, GridFunction &u, int depth_limit = 100);
+
+// Velocity coefficient
+void velocity_function(const Vector &x, Vector &v);
+
+// Initial condition
+double u0_function(const Vector &x, double);
+
+// Inflow boundary condition
+double inflow_function(const Vector &x);
+
+// Mesh bounding box
+Vector bb_min, bb_max;
+
+class FE_Evolution : public TimeDependentOperator
+{
+private:
+   BilinearForm &M, &K;
+   const Vector &b;
+   Solver *M_prec;
+   CGSolver M_solver;
+
+   mutable Vector z;
+
+public:
+   FE_Evolution(BilinearForm &M_, BilinearForm &K_, const Vector &b_);
+   void Update();
+   virtual void Mult(const Vector &x, Vector &y) const;
+
+   virtual ~FE_Evolution();
+};
+
+
+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 = 1;
+   sx = 1.0;
+   sy = 1.0;
+   double t_final = 1.0;
+   double dt = 0.002;
+   bool visualization = true;
+   int vis_steps = 5;
+
+   int precision = 8;
+   cout.precision(precision);
+
+   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.");
+   args.AddOption(&order, "-o", "--order",
+                  "Order (degree) of the finite elements.");
+   args.AddOption(&t_final, "-tf", "--t-final",
+                  "Final time; start time is 0.");
+   args.AddOption(&dt, "-dt", "--time-step",
+                  "Time step.");
+   args.AddOption(&sx, "-sx", "--sx",
+                  "mesh length in x direction");
+   args.AddOption(&sy, "-sy", "--sy",
+                  "mesh length in y direction");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&vis_steps, "-vs", "--visualization-steps",
+                  "Visualize every n-th timestep.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   Mesh mesh0 = Mesh::MakeCartesian2D(16, 16, mfem::Element::QUADRILATERAL,false,
+                                      sx,
+                                      sy);
+
+
+   std::vector<Vector> translations = {Vector({sx,0.0}), Vector({0.0,sy})};
+
+
+   Mesh mesh = Mesh::MakePeriodic(mesh0,
+                                  mesh0.CreatePeriodicVertexMapping(translations));
+
+
+   mesh.EnsureNCMesh();
+
+   int dim = mesh.Dimension();
+
+   ODESolver *ode_solver = new RK4Solver;
+
+
+   mesh.GetBoundingBox(bb_min, bb_max, max(order, 1));
+
+   // 5. Define the discontinuous DG finite element space of the given
+   //    polynomial order on the refined mesh.
+   DG_FECollection fec(order, dim, BasisType::GaussLobatto);
+   FiniteElementSpace fes(&mesh, &fec);
+   FiniteElementSpace fes_old(&mesh, &fec);
+
+   cout << "Number of unknowns: " << fes.GetVSize() << endl;
+
+   VectorFunctionCoefficient velocity(dim, velocity_function);
+   FunctionCoefficient inflow(inflow_function);
+   FunctionCoefficient u0(u0_function);
+
+   BilinearForm m(&fes);
+   BilinearForm k(&fes);
+   m.AddDomainIntegrator(new MassIntegrator);
+   constexpr double alpha = -1.0;
+   k.AddDomainIntegrator(new ConvectionIntegrator(velocity, alpha));
+   k.AddInteriorFaceIntegrator(
+      new NonconservativeDGTraceIntegrator(velocity, alpha));
+   k.AddBdrFaceIntegrator(
+      new NonconservativeDGTraceIntegrator(velocity, alpha));
+
+   LinearForm b(&fes);
+   b.AddBdrFaceIntegrator(
+      new BoundaryFlowIntegrator(inflow, velocity, alpha,-0.5));
+
+   m.Assemble();
+   int skip_zeros = 0;
+   k.Assemble(skip_zeros);
+   b.Assemble();
+   m.Finalize();
+   k.Finalize(skip_zeros);
+
+   // 7. Define the initial conditions, save the corresponding grid function to
+   //    a file and (optionally) save data in the VisIt format and initialize
+   //    GLVis visualization.
+   GridFunction u(&fes);
+   u0.SetTime(0.);
+   u.ProjectCoefficient(u0);
+
+
+   L2_FECollection orders_fec(0,dim);
+   FiniteElementSpace orders_fes(&mesh,&orders_fec);
+   GridFunction orders_gf(&orders_fes);
+   for (int i = 0; i<mesh.GetNE(); i++) { orders_gf(i) = order; }
+
+   socketstream sout;
+   // socketstream meshout;
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      sout.open(vishost, visport);
+      // meshout.open(vishost, visport);
+      if (!sout)
+      {
+         cout << "Unable to connect to GLVis server at "
+              << vishost << ':' << visport << endl;
+         visualization = false;
+         cout << "GLVis visualization disabled.\n";
+      }
+      else
+      {
+         sout.precision(precision);
+         sout << "solution\n" << mesh << u;
+         sout << flush;
+         cin.get();
+         // meshout.precision(precision);
+         // meshout << "solution\n" << mesh << orders_gf;
+         // meshout << "mesh\n" << mesh ;
+         // meshout << flush;
+         // cin.get();
+      }
+   }
+
+
+   // 8. Define the time-dependent evolution operator describing the ODE
+   //    right-hand side, and perform time-integration (looping over the time
+   //    iterations, ti, with a time-step dt).
+   FE_Evolution adv(m, k, b);
+
+   double t = 0.0;
+   adv.SetTime(t);
+
+   FunctionCoefficient u_ex(u0_function);
+
+   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;
+         u_ex.SetTime(t);
+         // Prefine(fes_old,u,u_ex, orders_gf, 5e-5, 5e-4);
+
+         mfem::out << "Global L2 Error = " << u.ComputeL2Error(u_ex) << std::endl;
+         Hrefine2(u,u_ex, 5e-5, 5e-4);
+         // mfem::out << "number of elements = " << mesh.GetNE() << endl;
+
+
+         m.Update();
+         m.Assemble();
+         m.Finalize();
+         k.Update();
+         k.Assemble(skip_zeros);
+         k.Finalize(skip_zeros);
+         b.Update();
+         b.Assemble();
+         adv.Update();
+         ode_solver->Init(adv);
+         if (visualization)
+         {
+            // GridFunction gf_ex(&fes);
+            // gf_ex.ProjectCoefficient(u_ex);
+            GridFunction * pr_u = ProlongToMaxOrder(&u);
+            sout << "solution\n" << mesh << *pr_u << flush;
+
+            // meshout << "solution\n" << mesh << orders_gf << flush;
+            // meshout << "mesh\n" << mesh << flush;
+         }
+      }
+   }
+
+   // 10. Free the used memory.
+   delete ode_solver;
+
+   return 0;
+}
+
+
+// Implementation of class FE_Evolution
+FE_Evolution::FE_Evolution(BilinearForm &M_, BilinearForm &K_, const Vector &b_)
+   : TimeDependentOperator(M_.Height()), M(M_), K(K_), b(b_), z(M_.Height())
+{
+   Array<int> ess_tdof_list;
+   M_prec = new OperatorJacobiSmoother(M, ess_tdof_list);
+   M_solver.SetPreconditioner(*M_prec);
+   M_solver.SetOperator(M);
+   M_solver.iterative_mode = false;
+   M_solver.SetRelTol(1e-9);
+   M_solver.SetAbsTol(0.0);
+   M_solver.SetMaxIter(100);
+   M_solver.SetPrintLevel(0);
+}
+
+void FE_Evolution::Update()
+{
+   height = M.Height();
+   width = M.Width();
+   z.SetSize(M.Height());
+
+   Array<int> ess_tdof_list;
+   delete M_prec;
+   M_prec = new OperatorJacobiSmoother(M, ess_tdof_list);
+   M_solver.SetPreconditioner(*M_prec);
+   M_solver.SetOperator(M);
+   M_solver.iterative_mode = false;
+   M_solver.SetRelTol(1e-9);
+   M_solver.SetAbsTol(0.0);
+   M_solver.SetMaxIter(100);
+   M_solver.SetPrintLevel(0);
+}
+
+void FE_Evolution::Mult(const Vector &x, Vector &y) const
+{
+   // y = M^{-1} (K x + b)
+   K.Mult(x, z);
+   z += b;
+   M_solver.Mult(z, y);
+}
+
+FE_Evolution::~FE_Evolution()
+{
+   delete M_prec;
+}
+
+// Velocity coefficient
+void velocity_function(const Vector &x, Vector &v)
+{
+   v.SetSize(2);
+   v(0) = 1.;
+   v(1) = 1.;
+}
+
+// Initial condition
+double u0_function(const Vector &x, double t)
+{
+   // give x0, y0;
+
+   // Rotation matrix
+   double theta = M_PI/4;
+   //
+
+   double x0 = 0.5;
+   double y0 = 0.5;
+   double w = 100.;
+   double c = 1.;
+   double ds = c*t;
+   Vector a(2);
+   a(0) = cos(theta);
+   a(1) = sin(theta);
+   // double xx = x(0) - a(0)*ds;
+   // double yy = x(1) - a(1)*ds;
+   double xx = x(0) - ds;
+   double yy = x(1) - ds;
+
+   double tol = 1e-6;
+   if (xx>= sx+tol || xx<= 0.0-tol)
+   {
+      xx -= floor(xx/sx) * sx;
+   }
+   if (yy>= sy+tol || yy<= 0.0-tol)
+   {
+      yy -= floor(yy/sy) * sy;
+   }
+
+
+
+   // double d = (xx-x0)*a(0) + (yy-y0)*a(1);
+   // double d1 = (xx-x0-0.5)*a(0) + (yy-y0-0.5)*a(1);
+   // double d2 = (xx-x0+0.5)*a(0) + (yy-y0+0.5)*a(1);
+   // return 1. + exp(-w*(d*d)) + exp(-w*(d1*d1)) + exp(-w*(d2*d2));
+   double dr_x = (xx-x0)*(xx-x0);
+   double dr_y = (yy-y0)*(yy-y0);
+   return 1. + exp(-w*(dr_x+dr_y));
+   // return 1. + exp(-w*(dr_x));
+}
+
+// Inflow boundary condition (zero for the problems considered in this example)
+double inflow_function(const Vector &x)
+{
+   return 1.0;
+}
+
+
+void Prefine(FiniteElementSpace & fes_old,
+             GridFunction &u, Coefficient &ex, GridFunction &orders_gf,
+             double min_thresh, double max_thresh)
+{
+   // get element errors
+   FiniteElementSpace * fes = u.FESpace();
+   int ne = fes->GetMesh()->GetNE();
+   Vector errors(ne);
+   u.ComputeElementL2Errors(ex,errors);
+   for (int i = 0; i<ne; i++)
+   {
+      double error = errors(i);
+      int order = fes->GetElementOrder(i);
+      if (error < min_thresh && order > 1)
+      {
+         fes->SetElementOrder(i,order-1);
+      }
+      else if (error > max_thresh && order < 2)
+      {
+         fes->SetElementOrder(i, order+1);
+      }
+      else
+      {
+         // do nothing
+      }
+   }
+
+   fes->Update(false);
+
+   PRefinementTransferOperator * T = new PRefinementTransferOperator(fes_old,*fes);
+
+   GridFunction u_fine(fes);
+   T->Mult(u,u_fine);
+
+   // copy the orders to the old space
+   for (int i = 0; i<ne; i++)
+   {
+      int order = fes->GetElementOrder(i);
+      fes_old.SetElementOrder(i,order);
+      orders_gf(i) = order;
+   }
+   fes_old.Update(false);
+
+   delete T;
+
+   // update old gridfuntion;
+   u = u_fine;
+
+}
+
+void Hrefine2(GridFunction &u, Coefficient & ex_coeff, double min_thresh,
+              double max_thresh)
+{
+   FiniteElementSpace * fes = u.FESpace();
+   Mesh * mesh = fes->GetMesh();
+   int ne = mesh->GetNE();
+   Vector errors(ne);
+   u.ComputeElementL2Errors(ex_coeff,errors);
+
+   Array<int> actions(ne);
+   for (int i = 0; i<ne; i++)
+   {
+      double error = errors(i);
+      if (error > max_thresh)
+      {
+         actions[i] = 1;
+      }
+      else if (error < min_thresh)
+      {
+         actions[i] = -1;
+      }
+      else
+      {
+         actions[i] = 0;
+      }
+   }
+   Refine(actions,u,1);
+
+   // construct a list of possible ref actions
+   // Array<int> actions(ne);
+   // for (int i = 0; i<ne; i++)
+   // {
+   //    double error = errors(i);
+   //    if (error > max_thresh && mesh->ncmesh->GetElementDepth(i) < 1)
+   //    {
+   //       actions[i] = 1;
+   //    }
+   //    else
+   //    {
+   //       actions[i] = 0;
+   //    }
+   // }
+
+   // // list of possible dref actions
+   // Array<int> derefactions(ne); derefactions = 0;
+   // const Table & dref_table = mesh->ncmesh->GetDerefinementTable();
+   // for (int i = 0; i<dref_table.Size(); i++)
+   // {
+   //    int size = dref_table.RowSize(i);
+   //    const int * row = dref_table.GetRow(i);
+   //    double error = 0.;
+   //    for (int j = 0; j<size; j++)
+   //    {
+   //       error += errors[row[j]];
+   //    }
+   //    if (error < min_thresh)
+   //    {
+   //       for (int j = 0; j<size; j++)
+   //       {
+   //          actions[row[j]] += -1;
+   //       }
+   //    }
+   // }
+
+   // // now refine the elements that have score >0 and deref the elements that have score < 0
+   // Array<Refinement> elements_to_refine;
+   // for (int i = 0; i<ne; i++)
+   // {
+   //    if (actions[i] > 0)
+   //    {
+   //       elements_to_refine.Append(Refinement(i,0b01));
+   //    }
+   // }
+
+   // mesh->GeneralRefinement(elements_to_refine);
+   // fes->Update();
+   // u.Update();
+
+   // // map old actions to new mesh
+   // Array<int> new_actions(mesh->GetNE());
+   // if (mesh->GetLastOperation() == mesh->REFINE)
+   // {
+   //    const CoarseFineTransformations &tr = mesh->GetRefinementTransforms();
+   //    Table coarse2fine;
+   //    tr.MakeCoarseToFineTable(coarse2fine);
+   //    new_actions = 1;
+   //    for (int i = 0; i<coarse2fine.Size(); i++)
+   //    {
+   //       if (coarse2fine.RowSize(i) == 1)
+   //       {
+   //          int * el = coarse2fine.GetRow(i);
+   //          new_actions[el[0]] = actions[i];
+   //       }
+   //    }
+   // }
+   // else
+   // {
+   //    new_actions = actions;
+   // }
+
+   // // create a dummy error vector
+   // Vector new_errors(mesh->GetNE());
+   // new_errors = infinity();
+   // for (int i = 0; i< new_errors.Size(); i++)
+   // {
+   //    if (new_actions[i] < 0)
+   //    {
+   //       new_errors[i] = 0.;
+   //    }
+   // }
+
+   // // any threshold would do here
+   // mesh->DerefineByError(new_errors,min_thresh);
+
+   // fes->Update();
+   // u.Update();
+}
+
+
+Table * Refine(Array<int> ref_actions, GridFunction &u, int depth_limit)
+{
+   FiniteElementSpace * fes = u.FESpace();
+   Mesh * mesh = fes->GetMesh();
+   int ne = mesh->GetNE();
+
+
+   //  ovewrite to no action if an element is marked for refinement but it exceeds the depth limit
+   for (int i = 0; i<ne; i++)
+   {
+      int depth = mesh->ncmesh->GetElementDepth(i);
+      if (depth >= depth_limit && ref_actions[i] == 1)
+      {
+         ref_actions[i] = 0;
+      }
+   }
+
+   // current policy to map agent_actions to actions
+   // 1. All elements that are marked for refinement are to perform the refinement
+   // 2. All of the "siblings" (i) of a marked element for refinement are assigned action=max(0,agent_actions[i])
+   //   i.e., a) if the action is to be refined then they are refined
+   //         b) if the action is to be derefined or no action then they get no action
+   // 3. If among the "siblings" there is no refinement action then the group is marked
+   //    for derefinement if the majority (including a tie) of the siblings are marked for derefinement
+   //    otherwise they are marked for no action
+   //  h-refine:   action =  1
+   //  h-derefine: action = -1
+   //  do nothing: action =  0
+
+   Array<int> actions(ne);
+   Array<int> actions_marker(ne);
+   actions_marker = 0;
+
+   const Table & deref_table = mesh->ncmesh->GetDerefinementTable();
+
+   for (int i = 0; i<deref_table.Size(); i++)
+   {
+      int n = deref_table.RowSize(i);
+      const int * row = deref_table.GetRow(i);
+      int sum_of_actions = 0;
+      bool ref_flag = false;
+      for (int j = 0; j<n; j++)
+      {
+         int action = ref_actions[row[j]];
+         sum_of_actions+=action;
+         if (action == 1)
+         {
+            ref_flag = true;
+            break;
+         }
+      }
+      if (ref_flag)
+      {
+         for (int j = 0; j<n; j++)
+         {
+            actions[row[j]] = max(0,ref_actions[row[j]]);
+            actions_marker[row[j]] = 1;
+         }
+      }
+      else
+      {
+         bool dref_flag = (2*abs(sum_of_actions) >= n) ? true : false;
+         for (int j = 0; j<n; j++)
+         {
+            actions[row[j]] = (dref_flag) ? -1 : 0;
+            actions_marker[row[j]] = 1;
+         }
+      }
+   }
+
+   for (int i = 0; i<ne; i++)
+   {
+      if (actions_marker[i] != 1)
+      {
+         if (ref_actions[i] == -1)
+         {
+            actions[i] = 0;
+         }
+         else
+         {
+            actions[i] = ref_actions[i];
+         }
+      }
+   }
+
+   // now the actions array holds feasible actions of -1,0,1
+   Array<Refinement> refinements;
+   for (int i = 0; i<ne; i++)
+   {
+      if (actions[i] == 1) {refinements.Append(Refinement(i,0b11));}
+   }
+   if (refinements.Size())
+   {
+      mesh->GeneralRefinement(refinements);
+      fes->Update();
+      u.Update();
+      ne = mesh->GetNE();
+   }
+
+   Table * ref_table = nullptr;
+   Table * dref_table = nullptr;
+   // now the derefinements
+   Array<int> new_actions(ne);
+   if (refinements.Size())
+   {
+      new_actions = 1;
+      const CoarseFineTransformations & tr = mesh->GetRefinementTransforms();
+      ref_table = new Table();
+      tr.MakeCoarseToFineTable(*ref_table);
+      for (int i = 0; i<ref_table->Size(); i++)
+      {
+         int n = ref_table->RowSize(i);
+         if (n == 1)
+         {
+            int * row = ref_table->GetRow(i);
+            new_actions[row[0]] = actions[i];
+         }
+      }
+   }
+   else
+   {
+      new_actions = actions;
+   }
+
+   Vector dummy_errors(ne);
+   dummy_errors = 1.0;
+   for (int i = 0; i<ne; i++)
+   {
+      if (new_actions[i] < 0)
+      {
+         dummy_errors[i] = 0.;
+      }
+   }
+   mesh->DerefineByError(dummy_errors,0.5);
+
+   fes->Update();
+   u.Update();
+
+   if (mesh->GetNE() < ne)
+   {
+      const CoarseFineTransformations & tr =
+         mesh->ncmesh->GetDerefinementTransforms();
+      Table coarse_to_fine_table;
+      tr.MakeCoarseToFineTable(coarse_to_fine_table);
+      dref_table = Transpose(coarse_to_fine_table);
+   }
+
+   // Build combined table of mesh modifications
+   Table * T = nullptr;
+   if (ref_table && dref_table)
+   {
+      T = Mult(*ref_table, * dref_table);
+      delete dref_table;
+      delete ref_table;
+   }
+   else if (ref_table)
+   {
+      T = ref_table;
+      delete dref_table;
+   }
+   else if (dref_table)
+   {
+      T= dref_table;
+      delete ref_table;
+   }
+   else
+   {
+      // do nothing: no mesh modifications happened
+   }
+
+   return T;
+}

examples/advection_href_reference_estimator.cpp

@@ -0,0 +1,957 @@
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+#include <algorithm>
+
+using namespace std;
+using namespace mfem;
+
+double sx, sy;
+
+Vector vel;
+
+
+void Prefine(FiniteElementSpace & fes_old,
+             GridFunction &u, GridFunction &pref_gf, GridFunction &orders_gf,
+             double min_thresh, double max_thresh);
+
+// void Prefine(FiniteElementSpace & fes_old,
+//              GridFunction &u, Coefficient &gf_ex, GridFunction &orders_gf,
+//              double min_thresh, double max_thresh);
+
+void Hrefine(GridFunction &u, Coefficient &gf_ex, double min_thresh,
+             double max_thresh);
+
+// void Hrefine2(GridFunction &u, Coefficient &gf_ex, double min_thresh,
+//   double max_thresh);
+Table * Hrefine2(GridFunction &u, GridFunction &u_ref, Table * refT,
+                 Coefficient &gf_ex, double min_thresh,
+                 double max_thresh);
+
+Table * Refine(Array<int> ref_actions, GridFunction &u, int depth_limit = 100);
+
+// Velocity coefficient
+void velocity_function(const Vector &x, Vector &v);
+
+// Initial condition
+double u0_function(const Vector &x, double);
+
+// Inflow boundary condition
+double inflow_function(const Vector &x);
+
+class FE_Evolution : public TimeDependentOperator
+{
+private:
+   BilinearForm &M, &K;
+   const Vector &b;
+   Solver *M_prec;
+   CGSolver M_solver;
+
+   mutable Vector z;
+
+public:
+   FE_Evolution(BilinearForm &M_, BilinearForm &K_, const Vector &b_);
+   void Update();
+   virtual void Mult(const Vector &x, Vector &y) const;
+
+   virtual ~FE_Evolution();
+};
+
+enum ref_kind
+{
+   order,        // p refinement
+   geometric     // h refinement
+};
+
+
+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 = 1;
+   sx = 1.0;
+   sy = 1.0;
+   double t_final = 1.0;
+   double dt = 0.002;
+   bool visualization = true;
+   int vis_steps = 5;
+   int refmode = 0;
+   int precision = 8;
+   ref_kind ref_mode = ref_kind::geometric;
+   cout.precision(precision);
+
+   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.");
+   args.AddOption(&order, "-o", "--order",
+                  "Order (degree) of the finite elements.");
+   args.AddOption(&t_final, "-tf", "--t-final",
+                  "Final time; start time is 0.");
+   args.AddOption(&dt, "-dt", "--time-step",
+                  "Time step.");
+   args.AddOption(&sx, "-sx", "--sx",
+                  "mesh length in x direction");
+   args.AddOption(&sy, "-sy", "--sy",
+                  "mesh length in y direction");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&vis_steps, "-vs", "--visualization-steps",
+                  "Visualize every n-th timestep.");
+   args.AddOption(&refmode, "-rm", "--refinement-mode",
+                  "0: 'p', 1: 'h' ");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   ref_mode = (ref_kind)refmode;
+
+   Mesh mesh0 = Mesh::MakeCartesian2D(32,32,mfem::Element::QUADRILATERAL,false,sx,
+                                      sy);
+   std::vector<Vector> translations = {Vector({sx,0.0}), Vector({0.0,sy})};
+   Mesh mesh = Mesh::MakePeriodic(mesh0,
+                                  mesh0.CreatePeriodicVertexMapping(translations));
+
+   mesh.EnsureNCMesh();
+
+   // compute reference solution
+   Mesh ref_mesh(mesh);
+
+   int dim = mesh.Dimension();
+
+   ODESolver *ode_solver = new RK4Solver;
+   ODESolver *ref_ode_solver = new RK4Solver;
+   ODESolver *pref_ode_solver = new RK4Solver;
+
+   // 5. Define the discontinuous DG finite element space of the given
+   //    polynomial order on the refined mesh.
+   DG_FECollection fec(order, dim, BasisType::GaussLobatto);
+   FiniteElementSpace fes(&mesh, &fec);
+   FiniteElementSpace ref_fes(&ref_mesh, &fec);
+   FiniteElementSpace fes_old(&mesh, &fec);
+   FiniteElementSpace pref_fes(&mesh, &fec);
+
+   for (int i = 0; i<mesh.GetNE(); i++)
+   {
+      int order = pref_fes.GetElementOrder(i);
+      pref_fes.SetElementOrder(i,order+2);
+   }
+   pref_fes.Update(false);
+
+
+   cout << "Number of unknowns: " << fes.GetVSize() << endl;
+
+   VectorFunctionCoefficient velocity(dim, velocity_function);
+   FunctionCoefficient inflow(inflow_function);
+   FunctionCoefficient u0(u0_function);
+
+   BilinearForm m(&fes);
+   BilinearForm k(&fes);
+   m.AddDomainIntegrator(new MassIntegrator);
+   constexpr double alpha = -1.0;
+   k.AddDomainIntegrator(new ConvectionIntegrator(velocity, alpha));
+   k.AddInteriorFaceIntegrator(
+      new NonconservativeDGTraceIntegrator(velocity, alpha));
+   k.AddBdrFaceIntegrator(
+      new NonconservativeDGTraceIntegrator(velocity, alpha));
+
+   LinearForm b(&fes);
+   b.AddBdrFaceIntegrator(
+      new BoundaryFlowIntegrator(inflow, velocity, alpha,-0.5));
+
+   m.Assemble();
+   int skip_zeros = 0;
+   k.Assemble(skip_zeros);
+   b.Assemble();
+   m.Finalize();
+   k.Finalize(skip_zeros);
+
+   // 7. Define the initial conditions, save the corresponding grid function to
+   //    a file and (optionally) save data in the VisIt format and initialize
+   //    GLVis visualization.
+   GridFunction u(&fes);
+   u0.SetTime(0.);
+   u.ProjectCoefficient(u0);
+
+   // reference solution (href)
+   BilinearForm m_ref(&ref_fes);
+   BilinearForm k_ref(&ref_fes);
+   m_ref.AddDomainIntegrator(new MassIntegrator);
+   k_ref.AddDomainIntegrator(new ConvectionIntegrator(velocity, alpha));
+   k_ref.AddInteriorFaceIntegrator(
+      new NonconservativeDGTraceIntegrator(velocity, alpha));
+   k_ref.AddBdrFaceIntegrator(
+      new NonconservativeDGTraceIntegrator(velocity, alpha));
+
+   LinearForm b_ref(&ref_fes);
+   b_ref.AddBdrFaceIntegrator(
+      new BoundaryFlowIntegrator(inflow, velocity, alpha,-0.5));
+
+   m_ref.Assemble();
+   k_ref.Assemble(skip_zeros);
+   b_ref.Assemble();
+   m_ref.Finalize();
+   k_ref.Finalize(skip_zeros);
+
+   // 7. Define the initial conditions, save the corresponding grid function to
+   //    a file and (optionally) save data in the VisIt format and initialize
+   //    GLVis visualization.
+   GridFunction u_ref(&ref_fes);
+   u_ref.ProjectCoefficient(u0);
+
+
+   Array<int> refinements(ref_mesh.GetNE());
+   refinements = 1;
+
+   Table * T1 = Refine(refinements, u_ref, 2);
+   refinements.SetSize(ref_mesh.GetNE());
+   refinements = 1;
+   Table * T2 = Refine(refinements, u_ref, 2);
+
+   Table * refT = Mult(*T1,*T2);
+
+   m_ref.Update();
+   m_ref.Assemble();
+   m_ref.Finalize();
+   k_ref.Update();
+   k_ref.Assemble(skip_zeros);
+   k_ref.Finalize(skip_zeros);
+   b_ref.Update();
+   b_ref.Assemble();
+
+
+
+   // reference solution (pref)
+   BilinearForm m_pref(&pref_fes);
+   BilinearForm k_pref(&pref_fes);
+   m_pref.AddDomainIntegrator(new MassIntegrator);
+   k_pref.AddDomainIntegrator(new ConvectionIntegrator(velocity, alpha));
+   k_pref.AddInteriorFaceIntegrator(
+      new NonconservativeDGTraceIntegrator(velocity, alpha));
+   k_pref.AddBdrFaceIntegrator(
+      new NonconservativeDGTraceIntegrator(velocity, alpha));
+
+   LinearForm b_pref(&pref_fes);
+   b_pref.AddBdrFaceIntegrator(
+      new BoundaryFlowIntegrator(inflow, velocity, alpha,-0.5));
+
+   m_pref.Assemble();
+   k_pref.Assemble(skip_zeros);
+   b_pref.Assemble();
+   m_pref.Finalize();
+   k_pref.Finalize(skip_zeros);
+
+   // 7. Define the initial conditions, save the corresponding grid function to
+   //    a file and (optionally) save data in the VisIt format and initialize
+   //    GLVis visualization.
+   GridFunction u_pref(&pref_fes);
+   u_pref.ProjectCoefficient(u0);
+   m_ref.Assemble();
+   m_ref.Finalize();
+   k_ref.Assemble(skip_zeros);
+   k_ref.Finalize(skip_zeros);
+
+
+   L2_FECollection orders_fec(0,dim);
+   FiniteElementSpace orders_fes(&mesh,&orders_fec);
+   GridFunction orders_gf(&orders_fes);
+   for (int i = 0; i<mesh.GetNE(); i++) { orders_gf(i) = order; }
+
+
+   socketstream sout;
+   // socketstream s_refout;
+   socketstream meshout;
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      sout.open(vishost, visport);
+      // s_refout.open(vishost, visport);
+      meshout.open(vishost, visport);
+      if (!sout)
+      {
+         cout << "Unable to connect to GLVis server at "
+              << vishost << ':' << visport << endl;
+         visualization = false;
+         cout << "GLVis visualization disabled.\n";
+      }
+      else
+      {
+         sout.precision(precision);
+         sout << "solution\n" << mesh << u;
+         sout << flush;
+         // s_refout.precision(precision);
+         // s_refout << "solution\n" << ref_mesh << u_ref;
+         // s_refout << flush;
+
+         meshout.precision(precision);
+         meshout << "solution\n" << mesh << orders_gf;
+         meshout << flush;
+         cin.get();
+      }
+   }
+
+
+   // 8. Define the time-dependent evolution operator describing the ODE
+   //    right-hand side, and perform time-integration (looping over the time
+   //    iterations, ti, with a time-step dt).
+   FE_Evolution adv(m, k, b);
+   FE_Evolution ref_adv(m_ref, k_ref, b_ref);
+   FE_Evolution pref_adv(m_pref, k_pref, b_pref);
+
+   double t = 0.0;
+   adv.SetTime(t);
+   double ref_t = 0.0;
+   ref_adv.SetTime(ref_t);
+
+   double pref_t = 0.0;
+   pref_adv.SetTime(pref_t);
+
+   FunctionCoefficient u_ex(u0_function);
+
+   ode_solver->Init(adv);
+   ref_ode_solver->Init(ref_adv);
+   pref_ode_solver->Init(pref_adv);
+
+   bool done = false;
+   for (int ti = 0; !done; )
+   {
+      double dt_real = min(dt, t_final - t);
+
+      ode_solver->Step(u, t, dt_real);
+      ref_ode_solver->Step(u_ref, ref_t, dt_real);
+      pref_ode_solver->Step(u_pref, pref_t, dt_real);
+      ti++;
+
+      done = (t >= t_final - 1e-8*dt);
+
+      if (done || ti % vis_steps == 0)
+      {
+         cout << "time step: " << ti << ", time: " << t << endl;
+         u_ex.SetTime(t);
+         // Prefine(fes_old,u,u_ex, orders_gf, 5e-5, 5e-4);
+
+         mfem::out << "Global L2 Error = " << u.ComputeL2Error(u_ex) << std::endl;
+         // Prefine(fes_old,u,u_ex, orders_gf, 5e-5, 5e-4);
+         // Hrefine2(u,u_ex, 5e-5, 5e-4);
+         // mfem::out << "refT size = " << refT->Size() << " x " << refT->Width() << endl;
+         if (ref_mode == ref_kind::geometric)
+         {
+            refT = Hrefine2(u,u_ref, refT, u_ex, 5e-5, 5e-4);
+         }
+         else
+         {
+            // Prefine(fes_old,u,u_ex, orders_gf, 5e-5, 5e-4);
+            Prefine(fes_old,u,u_pref, orders_gf, 5e-5, 5e-4);
+         }
+         // mfem::out << "refT size = " << refT->Size() << " x " << refT->Width() << endl;
+         // mfem::out << "number of elements = " << mesh.GetNE() << endl;
+
+
+         m.Update();
+         m.Assemble();
+         m.Finalize();
+         k.Update();
+         k.Assemble(skip_zeros);
+         k.Finalize(skip_zeros);
+         b.Update();
+         b.Assemble();
+         adv.Update();
+         ode_solver->Init(adv);
+
+
+         // ref_adv.Update();
+         // ref_ode_solver->Init(ref_adv);
+         if (visualization)
+         {
+            // GridFunction gf_ex(&fes);
+            // gf_ex.ProjectCoefficient(u_ex);
+            GridFunction * pr_u = ProlongToMaxOrder(&u);
+            sout << "solution\n" << mesh << *pr_u << flush;
+            // s_refout << "solution\n" << ref_mesh << u_ref << flush;
+
+            meshout << "solution\n" << mesh << orders_gf << flush;
+         }
+      }
+   }
+
+   // 10. Free the used memory.
+   delete ode_solver;
+   delete ref_ode_solver;
+   delete pref_ode_solver;
+
+   return 0;
+}
+
+
+// Implementation of class FE_Evolution
+FE_Evolution::FE_Evolution(BilinearForm &M_, BilinearForm &K_, const Vector &b_)
+   : TimeDependentOperator(M_.Height()), M(M_), K(K_), b(b_), z(M_.Height())
+{
+   Array<int> ess_tdof_list;
+   M_prec = new OperatorJacobiSmoother(M, ess_tdof_list);
+   M_solver.SetPreconditioner(*M_prec);
+   M_solver.SetOperator(M);
+   M_solver.iterative_mode = false;
+   M_solver.SetRelTol(1e-9);
+   M_solver.SetAbsTol(0.0);
+   M_solver.SetMaxIter(100);
+   M_solver.SetPrintLevel(0);
+}
+
+void FE_Evolution::Update()
+{
+   height = M.Height();
+   width = M.Width();
+   z.SetSize(M.Height());
+
+   Array<int> ess_tdof_list;
+   delete M_prec;
+   M_prec = new OperatorJacobiSmoother(M, ess_tdof_list);
+   M_solver.SetPreconditioner(*M_prec);
+   M_solver.SetOperator(M);
+   M_solver.iterative_mode = false;
+   M_solver.SetRelTol(1e-9);
+   M_solver.SetAbsTol(0.0);
+   M_solver.SetMaxIter(100);
+   M_solver.SetPrintLevel(0);
+}
+
+void FE_Evolution::Mult(const Vector &x, Vector &y) const
+{
+   // y = M^{-1} (K x + b)
+   K.Mult(x, z);
+   z += b;
+   M_solver.Mult(z, y);
+}
+
+FE_Evolution::~FE_Evolution()
+{
+   delete M_prec;
+}
+
+// Velocity coefficient
+void velocity_function(const Vector &x, Vector &v)
+{
+   v.SetSize(2);
+   v(0) = 1.;
+   v(1) = 1.;
+}
+
+// Initial condition
+double u0_function(const Vector &x, double t)
+{
+   // give x0, y0;
+
+   // Rotation matrix
+   double theta = M_PI/4;
+   //
+
+   double x0 = 0.5;
+   double y0 = 0.5;
+   double w = 100.;
+   double c = 1.;
+   double ds = c*t;
+   Vector a(2);
+   a(0) = cos(theta);
+   a(1) = sin(theta);
+   // double xx = x(0) - a(0)*ds;
+   // double yy = x(1) - a(1)*ds;
+   double xx = x(0) - ds;
+   double yy = x(1) - ds;
+
+   double tol = 1e-6;
+   if (xx>= sx+tol || xx<= 0.0-tol)
+   {
+      xx -= floor(xx/sx) * sx;
+   }
+   if (yy>= sy+tol || yy<= 0.0-tol)
+   {
+      yy -= floor(yy/sy) * sy;
+   }
+
+
+
+   // double d = (xx-x0)*a(0) + (yy-y0)*a(1);
+   // double d1 = (xx-x0-0.5)*a(0) + (yy-y0-0.5)*a(1);
+   // double d2 = (xx-x0+0.5)*a(0) + (yy-y0+0.5)*a(1);
+   // return 1. + exp(-w*(d*d)) + exp(-w*(d1*d1)) + exp(-w*(d2*d2));
+   double dr_x = (xx-x0)*(xx-x0);
+   double dr_y = (yy-y0)*(yy-y0);
+   return 1. + exp(-w*(dr_x+dr_y));
+   // return 1. + exp(-w*(dr_x));
+}
+
+// Inflow boundary condition (zero for the problems considered in this example)
+double inflow_function(const Vector &x)
+{
+   return 1.0;
+}
+
+
+// void Prefine(FiniteElementSpace & fes_old,
+//              GridFunction &u, Coefficient &ex, GridFunction &orders_gf,
+//              double min_thresh, double max_thresh)
+
+void Prefine(FiniteElementSpace & fes_old,
+             GridFunction &u, GridFunction &pref_gf, GridFunction &orders_gf,
+             double min_thresh, double max_thresh)
+{
+   // get element errors
+   FiniteElementSpace * fes = u.FESpace();
+   int ne = fes->GetMesh()->GetNE();
+   Vector errors(ne);
+
+   GridFunction pru(pref_gf.FESpace());
+   PRefinementTransferOperator * P = new PRefinementTransferOperator(*fes,
+                                                                     *pref_gf.FESpace());
+   P->Mult(u,pru);
+   delete P;
+   pru-=pref_gf;
+   ConstantCoefficient zero(0.0);
+   pru.ComputeElementL2Errors(zero,errors);
+
+   // u.ComputeElementL2Errors(ex,errors);
+   for (int i = 0; i<ne; i++)
+   {
+      double error = errors(i);
+      int order = fes->GetElementOrder(i);
+      if (error < min_thresh && order > 1)
+      {
+         fes->SetElementOrder(i,order-1);
+      }
+      else if (error > max_thresh && order < 2)
+      {
+         fes->SetElementOrder(i, order+1);
+      }
+      else
+      {
+         // do nothing
+      }
+   }
+
+   fes->Update(false);
+
+   PRefinementTransferOperator * T = new PRefinementTransferOperator(fes_old,*fes);
+
+   GridFunction u_fine(fes);
+   T->Mult(u,u_fine);
+
+   // copy the orders to the old space
+   for (int i = 0; i<ne; i++)
+   {
+      int order = fes->GetElementOrder(i);
+      fes_old.SetElementOrder(i,order);
+      orders_gf(i) = order;
+   }
+   fes_old.Update(false);
+
+   delete T;
+
+   // update old gridfuntion;
+   u = u_fine;
+
+}
+
+// void Hrefine2(GridFunction &u, Coefficient & ex_coeff, double min_thresh,
+//               double max_thresh)
+Table * Hrefine2(GridFunction &u, GridFunction &u_ref, Table * refT,
+                 Coefficient & ex_coeff, double min_thresh,
+                 double max_thresh)
+{
+   FiniteElementSpace * fes = u.FESpace();
+   Mesh * mesh = fes->GetMesh();
+   int ne = mesh->GetNE();
+   Vector errors(ne);
+   // u.ComputeElementL2Errors(ex_coeff,errors);
+
+   // copy the fespace, refine it up to element depth 2 and calculate the errors
+   // copy mesh
+   Mesh fine_mesh(*mesh);
+
+   FiniteElementSpace fes_copy(&fine_mesh,fes->FEColl());
+   GridFunction u_fine(&fes_copy);
+   // copy data;
+   u_fine = u;
+   Array<int>refinements(fine_mesh.GetNE());
+   refinements = 1;
+   Table * T1 = Refine(refinements,u_fine,2);
+   refinements.SetSize(fine_mesh.GetNE());
+   refinements = 1;
+   Table * T2 = Refine(refinements,u_fine,2);
+   Table * T = Mult(*T1, *T2);
+
+   delete T1;
+   delete T2;
+
+   // constract map
+   int n = T->Size();
+   int m = T->Width();
+   Array<int> elem_map(m);
+   for (int i = 0; i< n; i++)
+   {
+      int nr = T->RowSize(i);
+      int * row = T->GetRow(i);
+      int * ref_row = refT->GetRow(i);
+      for (int j = 0; j<nr ; j++ )
+      {
+         elem_map[row[j]] = ref_row[j];
+      }
+   }
+
+
+   // Table *fine2refT = Transpose(*Mult(*Transpose(*T), *refT));
+
+
+   // char vishost[] = "localhost";
+   // int  visport   = 19916;
+   // socketstream pr_out(vishost, visport);
+   // pr_out << "solution\n" << fine_mesh << u_fine << flush;
+
+   // calculate error
+
+   GridFunction diff(u_fine);
+   // this needs to change for reordering
+   diff-= u_ref;
+
+   ConstantCoefficient zero(0.0);
+   Vector fine_errors(fine_mesh.GetNE());
+   diff.ComputeElementL2Errors(zero,fine_errors);
+
+   // combine fine errors to current mesh;
+   // Table *Tt = Transpose(*T);
+   // mfem::out << "Tt->Size = " << Tt->Size() << endl;
+   // mfem::out << "errors = " << errors.Size() << endl;
+   for (int i = 0; i<T->Size(); i++)
+   {
+      int m = T->RowSize(i);
+      int *row = T->GetRow(i);
+      double err = 0.;
+      for (int j = 0; j<m; j++)
+      {
+         err += fine_errors[row[j]]*fine_errors[row[j]];
+      }
+      errors[i] = sqrt(err);
+   }
+
+
+
+   // cin.get();
+
+   // compute element errors by
+   // 1. refine the mesh up to mesh limit 2
+   // 2. Prolongate the current solution to the refined mesh
+   // 3. Calculate errors and combine them (for the coarse elements)
+   // 4. Derifine mesh
+
+   //copy the mesh
+   // Mesh * ref_mesh = new Mesh(*mesh);
+   // Array<int> ref_actions(ne);
+   // ref_actions = 1;
+   // NCMesh * ref_ncmesh = ref_mesh->ncmesh;
+
+
+
+
+   Array<int> actions(ne);
+   for (int i = 0; i<ne; i++)
+   {
+      double error = errors(i);
+      if (error > max_thresh)
+      {
+         actions[i] = 1;
+      }
+      else if (error < min_thresh)
+      {
+         actions[i] = -1;
+      }
+      else
+      {
+         actions[i] = 0;
+      }
+   }
+
+   Table * T3 = Refine(actions,u,1);
+
+
+   if (T3)
+   {
+      Table *Ttt = Mult(*Transpose(*T3), *refT);
+      delete T3;
+      delete refT;
+      refT = Ttt;
+   }
+
+   return refT;
+
+   // construct a list of possible ref actions
+   // Array<int> actions(ne);
+   // for (int i = 0; i<ne; i++)
+   // {
+   //    double error = errors(i);
+   //    if (error > max_thresh && mesh->ncmesh->GetElementDepth(i) < 1)
+   //    {
+   //       actions[i] = 1;
+   //    }
+   //    else
+   //    {
+   //       actions[i] = 0;
+   //    }
+   // }
+
+   // // list of possible dref actions
+   // Array<int> derefactions(ne); derefactions = 0;
+   // const Table & dref_table = mesh->ncmesh->GetDerefinementTable();
+   // for (int i = 0; i<dref_table.Size(); i++)
+   // {
+   //    int size = dref_table.RowSize(i);
+   //    const int * row = dref_table.GetRow(i);
+   //    double error = 0.;
+   //    for (int j = 0; j<size; j++)
+   //    {
+   //       error += errors[row[j]];
+   //    }
+   //    if (error < min_thresh)
+   //    {
+   //       for (int j = 0; j<size; j++)
+   //       {
+   //          actions[row[j]] += -1;
+   //       }
+   //    }
+   // }
+
+   // // now refine the elements that have score >0 and deref the elements that have score < 0
+   // Array<Refinement> elements_to_refine;
+   // for (int i = 0; i<ne; i++)
+   // {
+   //    if (actions[i] > 0)
+   //    {
+   //       elements_to_refine.Append(Refinement(i,0b01));
+   //    }
+   // }
+
+   // mesh->GeneralRefinement(elements_to_refine);
+   // fes->Update();
+   // u.Update();
+
+   // // map old actions to new mesh
+   // Array<int> new_actions(mesh->GetNE());
+   // if (mesh->GetLastOperation() == mesh->REFINE)
+   // {
+   //    const CoarseFineTransformations &tr = mesh->GetRefinementTransforms();
+   //    Table coarse2fine;
+   //    tr.MakeCoarseToFineTable(coarse2fine);
+   //    new_actions = 1;
+   //    for (int i = 0; i<coarse2fine.Size(); i++)
+   //    {
+   //       if (coarse2fine.RowSize(i) == 1)
+   //       {
+   //          int * el = coarse2fine.GetRow(i);
+   //          new_actions[el[0]] = actions[i];
+   //       }
+   //    }
+   // }
+   // else
+   // {
+   //    new_actions = actions;
+   // }
+
+   // // create a dummy error vector
+   // Vector new_errors(mesh->GetNE());
+   // new_errors = infinity();
+   // for (int i = 0; i< new_errors.Size(); i++)
+   // {
+   //    if (new_actions[i] < 0)
+   //    {
+   //       new_errors[i] = 0.;
+   //    }
+   // }
+
+   // // any threshold would do here
+   // mesh->DerefineByError(new_errors,min_thresh);
+
+   // fes->Update();
+   // u.Update();
+}
+
+
+Table * Refine(Array<int> ref_actions, GridFunction &u, int depth_limit)
+{
+   FiniteElementSpace * fes = u.FESpace();
+   Mesh * mesh = fes->GetMesh();
+   int ne = mesh->GetNE();
+
+
+   //  ovewrite to no action if an element is marked for refinement but it exceeds the depth limit
+   for (int i = 0; i<ne; i++)
+   {
+      int depth = mesh->ncmesh->GetElementDepth(i);
+      if (depth >= depth_limit && ref_actions[i] == 1)
+      {
+         ref_actions[i] = 0;
+      }
+   }
+
+   // current policy to map agent_actions to actions
+   // 1. All elements that are marked for refinement are to perform the refinement
+   // 2. All of the "siblings" (i) of a marked element for refinement are assigned action=max(0,agent_actions[i])
+   //   i.e., a) if the action is to be refined then they are refined
+   //         b) if the action is to be derefined or no action then they get no action
+   // 3. If among the "siblings" there is no refinement action then the group is marked
+   //    for derefinement if the majority (including a tie) of the siblings are marked for derefinement
+   //    otherwise they are marked for no action
+   //  h-refine:   action =  1
+   //  h-derefine: action = -1
+   //  do nothing: action =  0
+
+   Array<int> actions(ne);
+   Array<int> actions_marker(ne);
+   actions_marker = 0;
+
+   const Table & deref_table = mesh->ncmesh->GetDerefinementTable();
+
+   for (int i = 0; i<deref_table.Size(); i++)
+   {
+      int n = deref_table.RowSize(i);
+      const int * row = deref_table.GetRow(i);
+      int sum_of_actions = 0;
+      bool ref_flag = false;
+      for (int j = 0; j<n; j++)
+      {
+         int action = ref_actions[row[j]];
+         sum_of_actions+=action;
+         if (action == 1)
+         {
+            ref_flag = true;
+            break;
+         }
+      }
+      if (ref_flag)
+      {
+         for (int j = 0; j<n; j++)
+         {
+            actions[row[j]] = max(0,ref_actions[row[j]]);
+            actions_marker[row[j]] = 1;
+         }
+      }
+      else
+      {
+         bool dref_flag = (2*abs(sum_of_actions) >= n) ? true : false;
+         for (int j = 0; j<n; j++)
+         {
+            actions[row[j]] = (dref_flag) ? -1 : 0;
+            actions_marker[row[j]] = 1;
+         }
+      }
+   }
+
+   for (int i = 0; i<ne; i++)
+   {
+      if (actions_marker[i] != 1)
+      {
+         if (ref_actions[i] == -1)
+         {
+            actions[i] = 0;
+         }
+         else
+         {
+            actions[i] = ref_actions[i];
+         }
+      }
+   }
+
+   // now the actions array holds feasible actions of -1,0,1
+   Array<Refinement> refinements;
+   for (int i = 0; i<ne; i++)
+   {
+      if (actions[i] == 1) {refinements.Append(Refinement(i,0b11));}
+   }
+   if (refinements.Size())
+   {
+      mesh->GeneralRefinement(refinements);
+      fes->Update();
+      u.Update();
+      ne = mesh->GetNE();
+   }
+
+   Table * ref_table = nullptr;
+   Table * dref_table = nullptr;
+   // now the derefinements
+   Array<int> new_actions(ne);
+   if (refinements.Size())
+   {
+      new_actions = 1;
+      const CoarseFineTransformations & tr = mesh->GetRefinementTransforms();
+      ref_table = new Table();
+      tr.MakeCoarseToFineTable(*ref_table);
+      for (int i = 0; i<ref_table->Size(); i++)
+      {
+         int n = ref_table->RowSize(i);
+         if (n == 1)
+         {
+            int * row = ref_table->GetRow(i);
+            new_actions[row[0]] = actions[i];
+         }
+      }
+   }
+   else
+   {
+      new_actions = actions;
+   }
+
+   Vector dummy_errors(ne);
+   dummy_errors = 1.0;
+   for (int i = 0; i<ne; i++)
+   {
+      if (new_actions[i] < 0)
+      {
+         dummy_errors[i] = 0.;
+      }
+   }
+   mesh->DerefineByError(dummy_errors,0.5);
+
+   fes->Update();
+   u.Update();
+
+   if (mesh->GetNE() < ne)
+   {
+      const CoarseFineTransformations & tr =
+         mesh->ncmesh->GetDerefinementTransforms();
+      Table coarse_to_fine_table;
+      tr.MakeCoarseToFineTable(coarse_to_fine_table);
+      dref_table = Transpose(coarse_to_fine_table);
+   }
+
+   // Build combined table of mesh modifications
+   Table * T = nullptr;
+   if (ref_table && dref_table)
+   {
+      T = Mult(*ref_table, * dref_table);
+      delete dref_table;
+      delete ref_table;
+   }
+   else if (ref_table)
+   {
+      T = ref_table;
+      delete dref_table;
+   }
+   else if (dref_table)
+   {
+      T= dref_table;
+      delete ref_table;
+   }
+   else
+   {
+      // do nothing: no mesh modifications happened
+   }
+
+   return T;
+}

examples/ex18.cpp

@@ -45,7 +45,7 @@
 
 // Classes FE_Evolution, RiemannSolver, DomainIntegrator and FaceIntegrator
 // shared between the serial and parallel version of the example.
-#include "ex18.hpp"
+#include "fem/auxiliary.hpp"
 
 // Choice for the problem setup. See InitialCondition in ex18.hpp.
 int problem;
@@ -53,21 +53,84 @@ int problem;
 // Equation constant parameters.
 const int num_equation = 4;
 const double specific_heat_ratio = 1.4;
-const double gas_constant = 1.0;
+const double gas_constant = 8.3145;
 
 // Maximum characteristic speed (updated by integrators)
 double max_char_speed;
 
+// Initial condition
+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 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.
    problem = 1;
    const char *mesh_file = "../data/periodic-square.mesh";
    int ref_levels = 1;
-   int order = 3;
+   int order = 2;
    int ode_solver_type = 4;
    double t_final = 2.0;
-   double dt = -0.01;
+   double dt = 0.001;
    double cfl = 0.3;
    bool visualization = true;
    int vis_steps = 50;
@@ -189,17 +252,17 @@ 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 DomainIntegrator(dim));
+   Aflux.AddDomainIntegrator(new TransposeIntegrator(new GradientIntegrator()));
    Aflux.Assemble();
 
    NonlinearForm A(&vfes);
-   RiemannSolver rsolver;
-   A.AddInteriorFaceIntegrator(new FaceIntegrator(rsolver, dim));
+   RiemannSolver rsolver(specific_heat_ratio, num_equation);
+   A.AddInteriorFaceIntegrator(new FaceIntegrator(rsolver, dim, num_equation));
 
    // 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());
+   EulerSystem euler(vfes, A, Aflux.SpMat(), specific_heat_ratio, num_equation);
 
    // Visualize the density
    socketstream sout;
@@ -245,17 +308,6 @@ int main(int argc, char *argv[])
    double t = 0.0;
    euler.SetTime(t);
    ode_solver->Init(euler);
-
-   if (cfl > 0)
-   {
-      // Find a safe dt, using a temporary vector. Calling Mult() computes the
-      // maximum char speed at all quadrature points on all faces.
-      Vector z(A.Width());
-      max_char_speed = 0.;
-      A.Mult(sol, z);
-      dt = cfl * hmin / max_char_speed / (2*order+1);
-   }
-
    // Integrate in time.
    bool done = false;
    for (int ti = 0; !done; )
@@ -263,10 +315,6 @@ int main(int argc, char *argv[])
       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);
-      }
       ti++;
 
       done = (t >= t_final - 1e-8*dt);

examples/ref_tables.cpp

@@ -0,0 +1,220 @@
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+#include <algorithm>
+
+using namespace std;
+using namespace mfem;
+
+
+
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   bool visualization = true;
+
+   int precision = 8;
+   cout.precision(precision);
+
+   OptionsParser args(argc, argv);
+   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);
+
+
+   Mesh mesh = Mesh::MakeCartesian2D(2,2,mfem::Element::QUADRILATERAL);
+   mesh.EnsureNCMesh();
+
+   Array<Table * > map_table;
+   // in case of a refinement we need the Transpose of the coarse2fine Table
+   // in case of a derefinement we need the coarse2fine table
+
+   int ref = 2;
+   char vishost[] = "localhost";
+   int  visport   = 19916;
+
+   if (visualization)
+   {
+      socketstream sout(vishost, visport);
+      sout.precision(precision);
+      sout << "mesh\n" << mesh << flush;
+   }
+
+   Table temp;
+
+   Array<int> refinements;
+   refinements.Append(0);
+   refinements.Append(1);
+   refinements.Append(3);
+   mesh.GeneralRefinement(refinements);
+   if (visualization)
+   {
+      socketstream sout(vishost, visport);
+      sout.precision(precision);
+      sout << "mesh\n" << mesh << flush;
+   }
+   const CoarseFineTransformations & tr = mesh.GetRefinementTransforms();
+   tr.MakeCoarseToFineTable(temp);
+   map_table.Append(Transpose(temp));
+
+   refinements.SetSize(0);
+   refinements.Append(7);
+   mesh.GeneralRefinement(refinements);
+   if (visualization)
+   {
+      socketstream sout(vishost, visport);
+      sout.precision(precision);
+      sout << "mesh\n" << mesh << flush;
+   }
+   const CoarseFineTransformations & tr1 = mesh.GetRefinementTransforms();
+   tr1.MakeCoarseToFineTable(temp);
+   map_table.Append(Transpose(temp));
+
+
+   // // for (int i = 0; i<ref; i++)
+   // // {
+   // //    mesh.RandomRefinement(0.5);
+   // //    if (visualization)
+   // //    {
+   // //       socketstream sout(vishost, visport);
+   // //       sout.precision(precision);
+   // //       sout << "mesh\n" << mesh << flush;
+   // //    }
+   // //    const CoarseFineTransformations & tr = mesh.GetRefinementTransforms();
+   // //    tr.MakeCoarseToFineTable(temp);
+   // //    map_table.Append(Transpose(temp));
+   // // }
+
+   // // derefine
+   Vector errors(mesh.GetNE());
+   errors = 1;
+   errors[12] = 0.;  errors[14] = 0.;
+   errors[13] = 0.;  errors[15] = 0.;
+   // errors[21] = 0.; errors[18] = 0.;
+   // errors[19] = 0.; errors[20] = 0.;
+   mesh.DerefineByError(errors,0.3);
+
+   if (visualization)
+   {
+      socketstream sout1(vishost, visport);
+      sout1.precision(precision);
+      sout1 << "mesh\n" << mesh << flush;
+   }
+
+   const CoarseFineTransformations &tr2 = mesh.ncmesh->GetDerefinementTransforms();
+   map_table.Append(new Table);
+   tr2.MakeCoarseToFineTable(*map_table.Last());
+
+   for (int j = 0; j<mesh.GetNE(); j++)
+   {
+      // while the element depth > 0 for the element
+      mfem::out << "Refinement history for element: " << std::setw(4) << j << ": ";
+      int row = j;
+      for (int i = map_table.Size()-1 ; i>=0; i--)
+      {
+         row = map_table[i]->GetRow(row)[0];
+         if (i == 0)
+         {
+            mfem::out << std::setw(4) << row << endl;
+         }
+         else
+         {
+            mfem::out << std::setw(4) << row ;
+         }
+      }
+   }
+
+   // example ...
+   // create the maps from mesh 1 to mesh 3 (after one ref and 1 dref)
+   // 1. specify the newly created elements (either from ref or from dref)
+   // 2. mark the elements that are deleted
+   // 3. provide the map from the unmodified elements to their new numbers
+
+   // for the above scenario
+   Table & ref_table = *Transpose(*map_table[1]);
+   Table & dref_table = *Transpose(*map_table[2]);
+
+
+   Table * T = Mult(ref_table,dref_table);
+
+
+   mfem::out << "ref_table = " << endl;
+   ref_table.Print();
+
+   mfem::out << "dref_table = " << endl;
+   dref_table.Print();
+
+   mfem::out << "combined_table = " << endl;
+   T->Print();
+
+
+   Table * Tt = Transpose(*T);
+
+   Array<int> old_elems_map(T->Size()); // -1 if is to be deleted
+   Array<int> new_elems;
+
+   // loop though the old elements
+   for (int i = 0; i<T->Size(); i++)
+   {
+      // check row size
+      int n = T->RowSize(i);
+      int * row = T->GetRow(i);
+      if (n == 1)
+      {
+         // check the size of the transpose row
+         int m = Tt->RowSize(row[0]);
+         if (m == 1)
+         {
+            // the element is left unchanged
+            mfem::out << "Element number = " << i << " is mapped to element number = " <<
+                      row[0] << endl;
+            old_elems_map[i] = row[0];
+         }
+         else
+         {
+            mfem::out << "Element number = " << i <<
+                      " is derefined (deleted). Create new element: " << row[0] << endl;
+            old_elems_map[i] = -1;
+            new_elems.Append(row[0]);
+         }
+      }
+      else
+      {
+         mfem::out << "Element number = " << i <<
+                   " is refined (deleted). Create new elements = " ;
+         old_elems_map[i] = -1;
+         for (int j = 0; j<n; j++)
+         {
+            mfem::out << row[j];
+            new_elems.Append(row[j]);
+            if (j == n-1)
+            {
+               mfem::out << endl;
+            }
+            else
+            {
+               mfem::out << ", ";
+            }
+         }
+      }
+   }
+
+   // new_elems.Sort();
+   new_elems.Unique();
+
+   mfem::out << "elements map = " ; old_elems_map.Print(cout,old_elems_map.Size());
+   mfem::out << "new_elements = " ; new_elems.Print(cout,new_elems.Size());
+
+
+   return 0;
+}
+

fem/CMakeLists.txt

@@ -128,6 +128,7 @@ set(SRCS
   tmop_amr.cpp
   gslib.cpp
   transfer.cpp
+  auxiliary.cpp
   )
 
 set(HDRS
@@ -206,6 +207,7 @@ set(HDRS
   tmop_amr.hpp
   gslib.hpp
   transfer.hpp
+  auxiliary.hpp
   )
 
 if (MFEM_USE_SIDRE)

fem/auxiliary.hpp

@@ -0,0 +1,426 @@
+//                  MFEM Example 18 - Serial/Parallel Shared Code
+
+#include "mfem.hpp"
+using namespace std;
+using namespace mfem;
+
+
+// Time-dependent operator for the right-hand side of the ODE representing the
+// DG weak form.
+class EulerSystem : public TimeDependentOperator
+{
+private:
+   const int dim;
+   int num_equation;
+   double specific_heat_ratio;
+
+   FiniteElementSpace * vfes;
+   Operator &A;
+   SparseMatrix &Aflux;
+   std::vector<DenseMatrix> Me_inv;
+
+   mutable Vector state;
+   mutable DenseMatrix f;
+   mutable DenseTensor flux;
+   mutable Vector z;
+
+   void GetFlux(const DenseMatrix &state_, DenseTensor &flux_) const;
+
+public:
+   EulerSystem(FiniteElementSpace &vfes_,
+               Operator &A_, SparseMatrix &Aflux_,
+               double specific_heat_ratio_, int num_equation_);
+
+   virtual void Mult(const Vector &x, Vector &y) const;
+
+   virtual ~EulerSystem() {}
+
+};
+
+// Implements a simple Rusanov flux
+class RiemannSolver
+{
+public:
+   Vector flux1;
+   Vector flux2;
+
+   int num_equation;
+   double specific_heat_ratio;
+   RiemannSolver(double specific_heat_ratio_, int num_equation_);
+   double Eval(const Vector &state1, const Vector &state2,
+               const Vector &nor, Vector &flux);
+};
+
+// Interior face term: <F.n(u),[w]>
+class FaceIntegrator : public NonlinearFormIntegrator
+{
+public:
+   RiemannSolver rsolver;
+   int num_equation;
+   Vector shape1;
+   Vector shape2;
+   Vector funval1;
+   Vector funval2;
+   Vector nor;
+   Vector fluxN;
+
+   FaceIntegrator(RiemannSolver &rsolver_, const int dim, double num_equation_);
+
+   virtual void AssembleFaceVector(const FiniteElement &el1,
+                                   const FiniteElement &el2,
+                                   FaceElementTransformations &Tr,
+                                   const Vector &elfun, Vector &elvect);
+};
+
+// Implementation of class FE_Evolution
+EulerSystem::EulerSystem(FiniteElementSpace &vfes_,
+                         Operator &A_, SparseMatrix &Aflux_,
+                         double specific_heat_ratio_,
+                         int num_equation_)
+   : TimeDependentOperator(A_.Height()),
+     dim(vfes_.GetFE(0)->GetDim()),
+     vfes(&vfes_),
+     specific_heat_ratio(specific_heat_ratio_),
+     num_equation(num_equation_),
+     A(A_),
+     Aflux(Aflux_),
+     state(num_equation),
+     f(num_equation, dim),
+     flux(vfes->GetNDofs(), dim, num_equation),
+     z(A.Height())
+{
+   MassIntegrator mi;
+
+   for (int i = 0; i < vfes->GetNE(); i++) {
+      // Standard local assembly and inversion for energy mass matrices.
+      int dof = vfes->GetFE(i)->GetDof();
+      DenseMatrix Me(dof);
+      DenseMatrixInverse inv(&Me);
+
+      DenseMatrix inv_mi = DenseMatrix(vfes->GetFE(i)->GetDof(), vfes->GetFE(i)->GetDof());
+
+      mi.AssembleElementMatrix(*vfes->GetFE(i), *vfes->GetElementTransformation(i), Me);
+      inv.Factor();
+      inv.GetInverseMatrix(inv_mi);
+
+      Me_inv.push_back(inv_mi);
+   }
+}
+
+void EulerSystem::Mult(const Vector &x, Vector &y) const
+{
+   // 1. Create the vector z with the face terms -<F.n(u), [w]>.
+   A.Mult(x, z);
+
+   // 2. Add the element terms.
+   // i.  computing the flux approximately as a grid function by interpolating
+   //     at the solution nodes.
+   // ii. multiplying this grid function by a (constant) mixed bilinear form for
+   //     each of the num_equation, computing (F(u), grad(w)) for each equation.
+
+   DenseMatrix xmat(x.GetData(), vfes->GetNDofs(), num_equation);
+   GetFlux(xmat, flux);
+
+   for (int k = 0; k < num_equation; k++) {
+      Vector fk(flux(k).GetData(), dim * vfes->GetNDofs());
+      Vector zk(z.GetData() + k * vfes->GetNDofs(), vfes->GetNDofs());
+      Aflux.AddMult(fk, zk);
+   }
+
+   // 3. Multiply element-wise by the inverse mass matrices.
+   Vector zval;
+   Array<int> vdofs;
+
+   for (int i = 0; i < vfes->GetNE(); i++) {
+      int dof = vfes->GetFE(i)->GetDof();
+      DenseMatrix zmat, ymat(dof, num_equation);
+
+      // Return the vdofs ordered byNODES
+      vfes->GetElementVDofs(i, vdofs);
+      z.GetSubVector(vdofs, zval);
+      zmat.UseExternalData(zval.GetData(), dof, num_equation);
+      mfem::Mult(Me_inv[i], zmat, ymat);
+      y.SetSubVector(vdofs, ymat.GetData());
+   }
+}
+
+// Physicality check (at end)
+bool StateIsPhysical(const Vector &state, const int dim);
+
+// Pressure (EOS) computation
+inline double ComputePressure(const Vector &state, int num_equation,
+                              double specific_heat_ratio)
+{
+   const int udim = num_equation - 2;
+   const double den = state(0);
+   const Vector den_vel(state.GetData() + 1, udim);
+   const double den_energy = state(num_equation - 1);
+
+   double den_vel2 = 0;
+   for (int d = 0; d < udim; d++) {
+      den_vel2 += den_vel(d)*den_vel(d);
+   }
+   den_vel2 /= den;
+
+   return (specific_heat_ratio-1.0)*(den_energy - 0.5*den_vel2);
+}
+
+// Compute the vector flux F(u)
+void ComputeFlux(const Vector &state, int dim, DenseMatrix &flux, 
+                 double specific_heat_ratio, int num_equation)
+{
+   const int udim = num_equation - 2;
+   const double den = state(0);
+   const Vector den_vel(state.GetData() + 1, udim);
+   const double den_energy = state(num_equation - 1);
+
+   MFEM_ASSERT(StateIsPhysical(state, dim), "");
+
+   const double pres = ComputePressure(state, num_equation, specific_heat_ratio);
+   const double H = (den_energy + pres)/den;
+
+   // Hard-code quasi-1D cases
+   if (num_equation == 3) {
+      // Set x-flux
+      flux(0, 0) = den_vel(0);
+      flux(1, 0) = den_vel(0)*den_vel(0)/den + pres;
+      flux(2, 0) = den_vel(0)*H;
+
+      // Zero other components
+      for (int d = 1; d < dim; d++) {
+         for (int eq = 0; eq < num_equation; eq++){
+            flux(eq, d) = 0.0;
+         }
+      }
+   }
+   else {
+      MFEM_ASSERT(num_equation == dim + 2, "2D/3D solutions must be of size dim+2.")
+      for (int d = 0; d < dim; d++) {
+         flux(0, d) = den_vel(d);
+         for (int i = 0; i < dim; i++) {
+            flux(1+i, d) = den_vel(i) * den_vel(d) / den;
+         }
+         flux(1+d, d) += pres;
+      }
+      for (int d = 0; d < dim; d++) {
+         flux(1+dim, d) = den_vel(d) * H;
+      }
+   }
+}
+
+// Compute the scalar F(u).n
+void ComputeFluxDotN(const Vector &state, const Vector &nor, Vector &fluxN, 
+                     double specific_heat_ratio, int num_equation)
+{
+   const int udim = num_equation - 2;
+   // NOTE: nor in general is not a unit normal
+   const int dim = nor.Size();
+   MFEM_ASSERT(StateIsPhysical(state, dim), "");
+
+   DenseMatrix flux = DenseMatrix(num_equation, dim);
+   ComputeFlux(state, dim, flux, specific_heat_ratio, num_equation);
+
+   for (int i = 0; i < num_equation; i++) {
+      fluxN(i) = 0.0;
+      for (int d = 0; d < dim; d++) {
+         fluxN(i) += nor(d)*flux(i,d);
+      }
+   }
+}
+
+// Compute the maximum characteristic speed.
+inline double ComputeMaxCharSpeed(const Vector &state, const int dim,
+                                  double specific_heat_ratio, int num_equation)
+{
+   const int udim = num_equation - 2;
+   const double den = state(0);
+   const Vector den_vel(state.GetData() + 1, udim);
+
+   double den_vel2 = 0;
+   for (int d = 0; d < udim; d++) {
+      den_vel2 += den_vel(d)*den_vel(d);
+   }
+   den_vel2 /= den;
+
+   const double pres = ComputePressure(state, num_equation, specific_heat_ratio);
+   const double sound = sqrt(specific_heat_ratio*pres/den);
+   const double vel = sqrt(den_vel2/den);
+
+   return vel + sound;
+}
+
+// Compute the flux at solution nodes.
+void EulerSystem::GetFlux(const DenseMatrix &x_, DenseTensor &flux_) const
+{
+   const int flux_dof = flux_.SizeI();
+   const int flux_dim = flux_.SizeJ();
+
+   for (int i = 0; i < flux_dof; i++) {
+      for (int k = 0; k < num_equation; k++) {
+         state(k) = x_(i, k);
+      }
+
+      ComputeFlux(state, flux_dim, f, specific_heat_ratio, num_equation);
+
+      for (int d = 0; d < flux_dim; d++) {
+         for (int k = 0; k < num_equation; k++) {
+            flux_(i, d, k) = f(k, d);
+         }
+      }
+   }
+}
+
+// Implementation of class RiemannSolver
+RiemannSolver::RiemannSolver(double specific_heat_ratio_, int num_equation_) :
+   specific_heat_ratio(specific_heat_ratio_),
+   num_equation(num_equation_),
+   flux1(num_equation_),
+   flux2(num_equation_) { }
+
+double RiemannSolver::Eval(const Vector &state1, const Vector &state2,
+                           const Vector &nor, Vector &flux)
+{
+   // NOTE: nor in general is not a unit normal
+   const int dim = nor.Size();
+
+   MFEM_ASSERT(StateIsPhysical(state1, dim), "");
+   MFEM_ASSERT(StateIsPhysical(state2, dim), "");
+
+   const double maxE1 = ComputeMaxCharSpeed(state1, dim, specific_heat_ratio, num_equation);
+   const double maxE2 = ComputeMaxCharSpeed(state2, dim, specific_heat_ratio, num_equation);
+
+   const double maxE = max(maxE1, maxE2);
+
+   ComputeFluxDotN(state1, nor, flux1, specific_heat_ratio, num_equation);
+   ComputeFluxDotN(state2, nor, flux2, specific_heat_ratio, num_equation);
+
+   double normag = 0;
+   for (int i = 0; i < dim; i++) {
+      normag += nor(i) * nor(i);
+   }
+   normag = sqrt(normag);
+
+   for (int i = 0; i < num_equation; i++) {
+      flux(i) = 0.5 * (flux1(i) + flux2(i))
+                - 0.5 * maxE * (state2(i) - state1(i)) * normag;
+   }
+
+   return maxE;
+}
+
+
+// Implementation of class FaceIntegrator
+FaceIntegrator::FaceIntegrator(RiemannSolver &rsolver_, const int dim, double num_equation_) :
+   rsolver(rsolver_),
+   num_equation(num_equation_),
+   funval1(num_equation),
+   funval2(num_equation),
+   nor(dim),
+   fluxN(num_equation) { }
+
+void FaceIntegrator::AssembleFaceVector(const FiniteElement &el1,
+                                        const FiniteElement &el2,
+                                        FaceElementTransformations &Tr,
+                                        const Vector &elfun, Vector &elvect)
+{
+   // Compute the term <F.n(u),[w]> on the interior faces.
+   const int dof1 = el1.GetDof();
+   const int dof2 = el2.GetDof();
+
+   shape1.SetSize(dof1);
+   shape2.SetSize(dof2);
+
+   elvect.SetSize((dof1 + dof2) * num_equation);
+   elvect = 0.0;
+
+   DenseMatrix elfun1_mat(elfun.GetData(), dof1, num_equation);
+   DenseMatrix elfun2_mat(elfun.GetData() + dof1 * num_equation, dof2,
+                          num_equation);
+
+   DenseMatrix elvect1_mat(elvect.GetData(), dof1, num_equation);
+   DenseMatrix elvect2_mat(elvect.GetData() + dof1 * num_equation, dof2,
+                           num_equation);
+
+   // Integration order calculation from DGTraceIntegrator
+   int intorder;
+   if (Tr.Elem2No >= 0)
+      intorder = (min(Tr.Elem1->OrderW(), Tr.Elem2->OrderW()) +
+                  2*max(el1.GetOrder(), el2.GetOrder()));
+   else {
+      intorder = Tr.Elem1->OrderW() + 2*el1.GetOrder();
+   }
+
+   if (el1.Space() == FunctionSpace::Pk) {
+      intorder++;
+   }
+
+   const IntegrationRule *ir = &IntRules.Get(Tr.GetGeometryType(), intorder);
+   for (int i = 0; i < ir->GetNPoints(); i++) {
+      const IntegrationPoint &ip = ir->IntPoint(i);
+
+      Tr.SetAllIntPoints(&ip); // set face and element int. points
+
+      // Calculate basis functions on both elements at the face
+      el1.CalcShape(Tr.GetElement1IntPoint(), shape1);
+      el2.CalcShape(Tr.GetElement2IntPoint(), shape2);
+
+      // Interpolate elfun at the point
+      elfun1_mat.MultTranspose(shape1, funval1);
+      elfun2_mat.MultTranspose(shape2, funval2);
+
+      // Get the normal vector and the flux on the face
+      CalcOrtho(Tr.Jacobian(), nor);
+      const double mcs = rsolver.Eval(funval1, funval2, nor, fluxN);
+
+      fluxN *= ip.weight;
+      for (int k = 0; k < num_equation; k++) {
+         for (int s = 0; s < dof1; s++) {
+            elvect1_mat(s, k) -= fluxN(k) * shape1(s);
+         }
+         for (int s = 0; s < dof2; s++) {
+            elvect2_mat(s, k) += fluxN(k) * shape2(s);
+         }
+      }
+   }
+}
+
+// Check that the state is physical - enabled in debug mode
+bool StateIsPhysical(const Vector &state, const int dim)
+{
+   const double den = state(0);
+   const Vector den_vel(state.GetData() + 1, dim);
+   const double den_energy = state(1 + dim);
+
+   if (den < 0) {
+      cout << "Negative density: ";
+      for (int i = 0; i < state.Size(); i++) {
+         cout << state(i) << " ";
+      }
+      cout << endl;
+      return false;
+   }
+   if (den_energy <= 0) {
+      cout << "Negative energy: ";
+      for (int i = 0; i < state.Size(); i++) {
+         cout << state(i) << " ";
+      }
+      cout << endl;
+      return false;
+   }
+
+   double den_vel2 = 0;
+   for (int i = 0; i < dim; i++) { den_vel2 += den_vel(i) * den_vel(i); }
+   den_vel2 /= den;
+
+   const double int_energy = den_energy - 0.5 * den_vel2;
+
+   if (int_energy <= 0) {
+      cout << "Negative internal energy: " << int_energy << ", state: ";
+      for (int i = 0; i < state.Size(); i++) {
+         cout << state(i) << " ";
+      }
+      cout << endl;
+      return false;
+   }
+   return true;
+}

fem/gridfunc.cpp

@@ -4715,4 +4715,55 @@ GridFunction *Extrude1DGridFunction(Mesh *mesh, Mesh *mesh2d,
    return sol2d;
 }
 
+
+
+GridFunction* ProlongToMaxOrder(const GridFunction *x)
+{
+   const FiniteElementSpace *fespace = x->FESpace();
+   Mesh *mesh = fespace->GetMesh();
+   const FiniteElementCollection *fec = fespace->FEColl();
+
+   // find the max order in the space
+   int max_order = 1;
+   for (int i = 0; i < mesh->GetNE(); i++)
+   {
+      max_order = std::max(fespace->GetElementOrder(i), max_order);
+   }
+
+   // create a visualization space of max order for all elements
+   FiniteElementCollection *l2fec =
+      new L2_FECollection(max_order, mesh->Dimension(), BasisType::GaussLobatto);
+   FiniteElementSpace *l2space = new FiniteElementSpace(mesh, l2fec);
+
+   IsoparametricTransformation T;
+   DenseMatrix I;
+
+   GridFunction *prolonged_x = new GridFunction(l2space);
+
+   // interpolate solution vector in the larger space
+   for (int i = 0; i < mesh->GetNE(); i++)
+   {
+      Geometry::Type geom = mesh->GetElementGeometry(i);
+      T.SetIdentityTransformation(geom);
+
+      Array<int> dofs;
+      fespace->GetElementDofs(i, dofs);
+      Vector elemvect, l2vect;
+      x->GetSubVector(dofs, elemvect);
+
+      const auto *fe = fec->GetFE(geom, fespace->GetElementOrder(i));
+      const auto *l2fe = l2fec->GetFE(geom, max_order);
+
+      l2fe->GetTransferMatrix(*fe, T, I);
+      l2space->GetElementDofs(i, dofs);
+      l2vect.SetSize(dofs.Size());
+
+      I.Mult(elemvect, l2vect);
+      prolonged_x->SetSubVector(dofs, l2vect);
+   }
+
+   prolonged_x->MakeOwner(l2fec);
+   return prolonged_x;
+}
+
 }

fem/gridfunc.hpp

@@ -1101,6 +1101,9 @@ inline void QuadratureFunction::GetElementValues(int idx,
    }
 }
 
+
+GridFunction* ProlongToMaxOrder(const GridFunction *x);
+
 } // namespace mfem
 
 #endif

mesh/mesh.cpp

@@ -6156,6 +6156,59 @@ void Mesh::GetElementFaces(int i, Array<int> &el_faces, Array<int> &ori) const
    }
 }
 
+int Mesh::GetFaceElementsAndFaces(int face, Array<int> & elems,
+                                  Array<int> & faces) const
+{
+   int type = -1; // -1: bdr; 0: conforming, 1: slave, 2: master
+   bool nonconforming_face = ncmesh && (faces_info[face].NCFace != -1);
+   if (nonconforming_face)
+   {
+      int nc_index = faces_info[face].NCFace;
+      const NCFaceInfo &nc_info = nc_faces_info[nc_index];
+      const mfem::NCMesh::NCList &nc_list = ncmesh->GetNCList(Dim-1);
+
+      if (!nc_info.Slave)
+      {
+         type = 2;
+         elems.Append(ncmesh->elements[nc_list.masters[nc_index].element].index);
+         faces.Append(nc_list.masters[nc_index].index);
+         int j_begin = nc_list.masters[nc_index].slaves_begin;
+         int j_end = nc_list.masters[nc_index].slaves_end;
+         for (int j = j_begin; j<j_end ; j++)
+         {
+            elems.Append(ncmesh->elements[nc_list.slaves[j].element].index);
+            faces.Append(nc_list.slaves[j].index);
+         }
+         return type;
+      }
+      else
+      {
+         type = 1;
+         faces.Append(face);
+         faces.Append(nc_faces_info[nc_index].MasterFace);
+         int el1, el2;
+         GetFaceElements(face, &el1, &el2);
+         elems.Append(el1);
+         elems.Append(el2);
+         return type;
+      }
+   }
+   else
+   {
+      int el1, el2;
+      GetFaceElements(face, &el1, &el2);
+      elems.Append(el1);
+      faces.Append(face);
+      if (el2 != -1)
+      {
+         type = 0;
+         elems.Append(el2);
+         faces.Append(face);
+      }
+      return type;
+   }
+}
+
 void Mesh::GetBdrElementFace(int i, int *f, int *o) const
 {
    const int *bv, *fv;

mesh/mesh.hpp

@@ -1428,6 +1428,9 @@ public:
    void GetFaceInfos (int Face, int *Inf1, int *Inf2) const;
    void GetFaceInfos (int Face, int *Inf1, int *Inf2, int *NCFace) const;
 
+   int GetFaceElementsAndFaces(int face, Array<int> & elems,
+                               Array<int> & faces) const;
+
    Geometry::Type GetFaceGeometryType(int Face) const;
    Element::Type  GetFaceElementType(int Face) const;