vis script

.gitignore

@@ -52,6 +52,13 @@ examples/ex1[0-9]p
 examples/ex2[0-9]
 examples/ex2[0-9]p
 
+examples/soln-*
+examples/reference-*
+examples/*.dat
+examples/*.gif
+examples/*.mp4
+examples/*.gf
+
 examples/refined.mesh
 examples/displaced.mesh
 examples/mesh.*

examples/combine-figures

@@ -0,0 +1,30 @@
+#!/bin/bash
+
+# Use imagemagick to combine the figures into a composite view
+
+for i in {0..19}; do
+    ii=`printf %04d $i`
+
+    # label reference density
+    convert reference-$ii.png -gravity south -annotate -50+50 "Reference Fine" ,reference-$ii.png
+
+    # label optimal density
+    convert opt-soln-$ii.png -gravity south -annotate -50+50 "Optimal AMR" ,opt-soln-$ii.png
+
+    # label greedy density
+    convert greedy-soln-$ii.png -gravity south -annotate -50+50 "Max Current AMR" ,greedy-soln-$ii.png
+
+    # stack optimal on top of error plot
+    convert ,opt-soln-$ii.png opt-err-$i.png -append ,opt-soln-werr-$ii.png
+
+    # stack greedy on top of error plot
+    convert ,greedy-soln-$ii.png greedy-err-$i.png -append ,greedy-soln-werr-$ii.png
+    
+    # lay the 3 side by side
+    convert ,reference-$ii.png ,opt-soln-werr-$ii.png ,greedy-soln-werr-$ii.png +append ,ref-sbs-soln-$ii.1.png
+
+    # Label the composite image with timestep
+    convert ,ref-sbs-soln-$ii.1.png -gravity northwest -annotate +20+20 "Timestep $i" ref-sbs-soln-$ii.png
+done
+
+rm ,*

examples/ex18.cpp

@@ -42,6 +42,7 @@
 #include <fstream>
 #include <sstream>
 #include <iostream>
+#include <set>
 
 // Classes FE_Evolution, RiemannSolver, DomainIntegrator and FaceIntegrator
 // shared between the serial and parallel version of the example.
@@ -50,14 +51,177 @@
 // Choice for the problem setup. See InitialCondition in ex18.hpp.
 int problem;
 
+// The dim and hence number of equations is taken from the mesh.
+int dim;
+int num_equation;
+
 // Equation constant parameters.
-const int num_equation = 4;
 const double specific_heat_ratio = 1.4;
 const double gas_constant = 1.0;
 
 // Maximum characteristic speed (updated by integrators)
 double max_char_speed;
 
+void amr_update(FiniteElementSpace& fes,
+                FiniteElementSpace& dfes,
+                FiniteElementSpace& vfes,
+                BlockVector& u_block,
+                GridFunction& rho,
+                GridFunction& rho_u,
+                GridFunction& rho_e,
+                GridFunction& sol)
+{
+   fes.Update();
+   dfes.Update();
+   vfes.Update();
+
+   rho.Update();
+   rho_u.Update();
+   rho_e.Update();
+
+   Array<int> offsets(num_equation + 1);
+   for (int k = 0; k <= num_equation; k++) {
+      offsets[k] = k * vfes.GetNDofs();
+   }
+   u_block.Update(offsets);
+
+   Vector& sub0 = u_block.GetBlock(0);
+   Vector& sub1 = u_block.GetBlock(1);
+   Vector& sub2 = u_block.GetBlock(2);
+
+   sub0 = rho;
+   sub1 = rho_u;
+   sub2 = rho_e;
+
+   sol.MakeRef(&vfes, u_block.GetData());
+
+   rho.MakeRef(&fes, u_block.GetData() + offsets[0]);
+   rho_u.MakeRef(&dfes, u_block.GetData() + offsets[1]);
+   rho_e.MakeRef(&fes, u_block.GetData() + offsets[2]);
+}
+
+double estimate_dt(double cfl, Mesh& mesh,
+                   GridFunction& sol, NonlinearForm& A, int order)
+{
+   // Determine the minimum element size.
+   double hmin = 0.0;
+   hmin = mesh.GetElementSize(0, 1);
+   for (int i = 1; i < mesh.GetNE(); i++)
+   {
+      hmin = min(mesh.GetElementSize(i, 1), hmin);
+   }
+
+   double dt = -0.01;
+   // 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);
+
+   return dt;
+}
+
+void show_maps(Mesh& mesh, vector<int>& coarse_map, vector<int>& fine_map)
+{
+   for (int i = 0; i < mesh.GetNE(); i++) {
+      int depth = mesh.ncmesh->GetElementDepth(i);
+      if (depth == 0) {
+         printf("| %2d  ",i);
+      }
+      else {
+         printf("|%2d|%2d",i,i+1);
+         i++;
+      }
+   }
+   printf("| mesh i\n");
+
+   for (size_t i = 0; i < coarse_map.size(); i++) {
+      printf("| %2d  ",coarse_map[i]);
+   }
+   printf("| coarse map\n");
+
+   for (size_t i = 0; i < coarse_map.size(); i++) {
+      printf("| %2d  ",fine_map[i]);
+   }
+   printf("| fine map\n");
+
+   for (size_t i = 0; i < coarse_map.size(); i++) {
+      printf("| %2lu  ",i);
+   }
+   printf("| ref idx\n");
+
+   // some validation checks
+   for (size_t i = 0; i < coarse_map.size(); i++) {
+      //printf("checking reference index %d\n",i);
+      if (coarse_map[i] > -1) assert(  fine_map[i] == -1);
+      if (fine_map[i]   > -1) assert(coarse_map[i] == -1);
+      if (coarse_map[i] > -1) assert(mesh.ncmesh->GetElementDepth(coarse_map[i]) == 0);
+      if (fine_map[i]   > -1) assert(mesh.ncmesh->GetElementDepth(fine_map[i])   == 1);
+      if (fine_map[i]   > -1) assert(mesh.ncmesh->GetElementDepth(fine_map[i]+1) == 1);
+   }
+}
+
+void compute_reference_errors(int ti,
+                              FiniteElementCollection& fec,
+                              Mesh& mesh,
+                              GridFunction& den,
+                              const vector<int>& coarse_map,
+                              const vector<int>& fine_map,
+                              vector<double>& errors)
+{
+   // Create a temp copy of the mesh for refinement to reference resolution
+   Mesh tmp_mesh(mesh);
+   FiniteElementSpace tmp_fes(&tmp_mesh, &fec);
+   GridFunction tmp_den(&tmp_fes);
+   tmp_den = den; // copy data
+
+   // refine mesh and soln to reference resolution
+   Array<int> refs;
+   for (int i = 0; i < tmp_mesh.GetNE(); i++) {
+      if (tmp_mesh.ncmesh->GetElementDepth(i) == 0) {
+         refs.Append(i);
+      }
+   }
+   tmp_mesh.GeneralRefinement(refs);
+   tmp_den.FESpace()->Update();
+   tmp_den.Update();
+
+   // read in reference mesh
+   ostringstream fn;
+   fn << "reference-" << ti << ".mesh";
+   Mesh ref_mesh(fn.str().c_str());
+
+   // read in reference density
+   ostringstream sol_name;
+   sol_name << "reference-" << ti << "-" << 0 << ".gf";
+   ifstream gf_ifs(sol_name.str());
+   GridFunction ref_den(&ref_mesh, gf_ifs);
+
+   // integrate L2 errors
+   GridFunctionCoefficient ref_den_coeff(&ref_den);
+   Vector err(ref_mesh.GetNE());
+   tmp_den.ComputeElementL2Errors(ref_den_coeff, err);
+
+   {
+      ostringstream sol_name;
+      sol_name << "refined-current-rho-" << ti << ".gf";
+      ofstream sol_ofs(sol_name.str().c_str());
+      sol_ofs << tmp_den;
+   }
+
+
+   // sum L2 errors to reference mesh. The reference errors are always
+   // on a fully refined mesh.
+   int k = 0;
+   for (size_t i = 0; i < coarse_map.size(); i++) {
+      double e0 = err[k];
+      double e1 = err[k+1];
+      errors[i] = sqrt(e0*e0+e1*e1);
+      k += 2;
+   }
+}
+
 int main(int argc, char *argv[])
 {
    // 1. Parse command-line options.
@@ -69,8 +233,15 @@ int main(int argc, char *argv[])
    double t_final = 2.0;
    double dt = -0.01;
    double cfl = 0.3;
-   bool visualization = true;
+   bool visualization = false;
    int vis_steps = 50;
+   Array<int> rseq;
+   int nseq = 0;
+   int regrid_period = 1;
+   int greedy_refine = 0;
+   int lookahead_refine = 0;
+   int output_cycle_soln = 0;
+   int output_cycle_errors = 0;
 
    int precision = 8;
    cout.precision(precision);
@@ -98,6 +269,18 @@ int main(int argc, char *argv[])
                   "Enable or disable GLVis visualization.");
    args.AddOption(&vis_steps, "-vs", "--visualization-steps",
                   "Visualize every n-th timestep.");
+   args.AddOption(&rseq, "-rseq", "--refine-sequence",
+                  "Element sequence to refine.");
+   args.AddOption(&regrid_period, "-rp", "--regrid-period",
+                  "Timesteps per regrid.");
+   args.AddOption(&greedy_refine, "-gr", "--greedy-refine",
+                  "Use greedy refinement strategy.");
+   args.AddOption(&lookahead_refine, "-lr", "--lookahead-refine",
+                  "Use lookahead refinement strategy.");
+   args.AddOption(&output_cycle_soln, "-cs", "--output-cycle-solutions",
+                  "Output per-cycle solutions.");
+   args.AddOption(&output_cycle_errors, "-ce", "--output-cycle-errors",
+                  "Output per-cycle errors.");
 
    args.Parse();
    if (!args.Good())
@@ -107,12 +290,15 @@ int main(int argc, char *argv[])
    }
    args.PrintOptions(cout);
 
+   printf("dt = %f\n",dt);
+
    // 2. Read the mesh from the given mesh file. This example requires a 2D
    //    periodic mesh, such as ../data/periodic-square.mesh.
    Mesh mesh(mesh_file, 1, 1);
+   mesh.EnsureNCMesh();
    const int dim = mesh.Dimension();
 
-   MFEM_ASSERT(dim == 2, "Need a two-dimensional mesh for the problem definition");
+   num_equation = 2+dim;
 
    // 3. Define the ODE solver used for time integration. Several explicit
    //    Runge-Kutta methods are available.
@@ -162,7 +348,9 @@ int main(int argc, char *argv[])
    BlockVector u_block(offsets);
 
    // Momentum grid function on dfes for visualization.
-   GridFunction mom(&dfes, u_block.GetData() + offsets[1]);
+   GridFunction rho(&fes, u_block.GetData() + offsets[0]);
+   GridFunction rho_u(&dfes, u_block.GetData() + offsets[1]);
+   GridFunction rho_e(&fes, u_block.GetData() + offsets[2]);
 
    // Initialize the state.
    VectorFunctionCoefficient u0(num_equation, InitialCondition);
@@ -184,6 +372,7 @@ int main(int argc, char *argv[])
          sol_ofs.precision(precision);
          sol_ofs << uk;
       }
+
    }
 
    // 7. Set up the nonlinear form corresponding to the DG discretization of the
@@ -199,7 +388,7 @@ int main(int argc, char *argv[])
    // 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());
+   FE_Evolution euler(vfes, A, Aflux);
 
    // Visualize the density
    socketstream sout;
@@ -219,7 +408,7 @@ int main(int argc, char *argv[])
       else
       {
          sout.precision(precision);
-         sout << "solution\n" << mesh << mom;
+         sout << "solution\n" << mesh << rho;
          sout << "pause\n";
          sout << flush;
          cout << "GLVis visualization paused."
@@ -227,17 +416,6 @@ int main(int argc, char *argv[])
       }
    }
 
-   // Determine the minimum element size.
-   double hmin = 0.0;
-   if (cfl > 0)
-   {
-      hmin = mesh.GetElementSize(0, 1);
-      for (int i = 1; i < mesh.GetNE(); i++)
-      {
-         hmin = min(mesh.GetElementSize(i, 1), hmin);
-      }
-   }
-
    // Start the timer.
    tic_toc.Clear();
    tic_toc.Start();
@@ -246,27 +424,348 @@ int main(int argc, char *argv[])
    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);
+   // map from base coarse numbering to current coarse numbering.
+   // -1 entries for elements that are currently refined.
+   vector<int> coarse_map(mesh.GetNE());
+   for (size_t i = 0; i < coarse_map.size(); ++i) {
+      coarse_map[i] = i;
    }
 
+   // map from base coarse numbering to current fine numbering. The
+   // second element is always this index+1.
+   // -1 entries for elements that are currently coarse.
+   vector<int> fine_map(mesh.GetNE());
+   for (size_t i = 0; i < fine_map.size(); ++i) {
+      fine_map[i] = -1;
+   }
+
+   //show_maps(mesh,coarse_map,fine_map);
+
+   // vector<int> coarse_map(mesh.GetNE());
+   // for (int i = 0; i < coarse_map.size(); ++i) {
+   //    coarse_map[i] = i;
+   // }
+
+   // the elements that are currently refined
+   set<int> cur_ref_set; // in base numbering
+
    // Integrate in time.
    bool done = false;
    for (int ti = 0; !done; )
    {
+
+      if (greedy_refine && !(ti % regrid_period )) {
+         vector<double> ref_errors(coarse_map.size());
+         compute_reference_errors(ti, fec, mesh, rho,
+                                  coarse_map, fine_map, ref_errors);
+         vector<double>::iterator it;
+         it = std::max_element(ref_errors.begin(), ref_errors.end());
+         int ref_max = std::distance(ref_errors.begin(), it);
+         std::cout << "max at: " << std::distance(ref_errors.begin(), it) << '\n';
+         rseq.SetSize(nseq+1);
+         rseq[nseq] = ref_max;
+      }
+
+      if (lookahead_refine && !(ti % regrid_period )) {
+
+         MFEM_ASSERT(cfl < 0, "need fixed time step for lookahead.");
+
+         double dt_real = min(dt, t_final - t);
+
+         // advance a copy of the solution
+         GridFunction tmp_sol(sol);
+
+         ode_solver->Step(tmp_sol, t, dt_real);
+
+         GridFunction tmp_rho(&fes, tmp_sol);
+
+         vector<double> ref_errors(coarse_map.size());
+         compute_reference_errors(ti+1, fec, mesh, tmp_rho,
+                                  coarse_map, fine_map, ref_errors);
+
+         vector<double>::iterator it;
+         it = std::max_element(ref_errors.begin(), ref_errors.end());
+         int ref_max = std::distance(ref_errors.begin(), it);
+         std::cout << "max at: " << std::distance(ref_errors.begin(), it) << '\n';
+         rseq.SetSize(nseq+1);
+         rseq[nseq] = ref_max;
+      }
+
+      // adapt mesh to new element(s)
+      if (rseq.Size() && !(ti % regrid_period )) {
+
+         //printf("*** begin refinement ***\n");
+         set<int> new_ref_set; // in base numbering
+
+         // use base numbering
+         int ne = coarse_map.size();
+         int el1 = rseq[nseq++];
+         int el2 = el1-1;
+         int el3 = el1+1;
+         el2 = (el2 + ne) % ne;
+         el3 = (el3 + ne) % ne;
+         new_ref_set.insert(el1);
+         new_ref_set.insert(el2);
+         new_ref_set.insert(el3);
+
+         // The coarsen set is elements from cur_ref_set not in new_ref_set
+         set<int> coarsen_set;
+         std::set_difference(cur_ref_set.begin(), cur_ref_set.end(),
+                             new_ref_set.begin(), new_ref_set.end(),
+                             std::inserter(coarsen_set, coarsen_set.begin()));
+
+         // The refine set is elements from new_ref_set not in cur_ref_set
+         set<int> refine_set;
+         std::set_difference(new_ref_set.begin(), new_ref_set.end(),
+                             cur_ref_set.begin(), cur_ref_set.end(),
+                             std::inserter(refine_set, refine_set.begin()));
+
+         set<int>::iterator it;
+
+         // printf("cur_ref set\n");
+         //
+         // for (it = cur_ref_set.begin(); it != cur_ref_set.end(); ++it) {
+         //    const int& i = *it;
+         //    printf("%d ",i);
+         // }
+         // printf("\n");
+
+         // printf("new_ref set\n");
+         // for (it = new_ref_set.begin(); it != new_ref_set.end(); ++it) {
+         //    const int& i = *it;
+         //    printf("%d ",i);
+         // }
+         // printf("\n");
+
+         // printf("refine set: ");
+         // for (it = refine_set.begin(); it != refine_set.end(); ++it) {
+         //    const int& i = *it;
+         //    printf("%d ",i);
+         // }
+         // printf("\n");
+
+         // printf("coarsen set: ");
+         // for (it = coarsen_set.begin(); it != coarsen_set.end(); ++it) {
+         //    const int& i = *it;
+         //    printf("%d ",i);
+         // }
+         // printf("\n");
+
+
+         // Perform any new refinements. Translate into current numbering.
+         Array<int> els;
+         for (it = refine_set.begin(); it != refine_set.end(); ++it) {
+            const int& i = *it;
+            assert(coarse_map[i] >= 0);
+            //printf("adding %d -> %d to be refined\n",i,coarse_map[i]);
+            els.Append(coarse_map[i]);
+         }
+
+         if (els.Size()) {
+            //printf("  ** start refinements **\n");
+            mesh.GeneralRefinement(els);
+
+            const CoarseFineTransformations &cft = mesh.GetRefinementTransforms();
+            //printf("updating coarse_map and fine_map after refinement\n");
+            Table c2f;
+            cft.GetCoarseToFineMap(mesh, c2f);
+            Array<int> row;
+            for (size_t i = 0; i < coarse_map.size(); ++i) {
+
+               int old = coarse_map[i];
+               if (old > -1) {
+                  // was coarse
+                  c2f.GetRow(old,row);
+                  if (row.Size() == 1) {
+                     // was coarse, still coarse
+                     coarse_map[i] = row[0];
+                  }
+                  else {
+                     // was coarse, now fine
+                     coarse_map[i] = -1;
+                     fine_map[i] = row[0];
+                  }
+               }
+               else {
+                  // was fine, still fine
+                  old = fine_map[i];
+                  assert(old != -1);
+                  c2f.GetRow(old,row);
+                  assert(row.Size() == 1);
+                  fine_map[i] = row[0];
+               }
+            }
+
+            // printf("new coarse_map is:\n");
+            // for (size_t i = 0; i < coarse_map.size(); i++) {
+            //    printf("%lu -> %d\n",i,coarse_map[i]);
+            // }
+            // printf("new fine_map is:\n");
+            // for (size_t i = 0; i < fine_map.size(); i++) {
+            //    printf("%lu -> %d\n",i,fine_map[i]);
+            // }
+
+            amr_update(fes, dfes, vfes, u_block, rho, rho_u, rho_e, sol);
+            A.Update();
+            Aflux.Update();
+            Aflux.Assemble();
+            euler.Update();
+            ode_solver->Init(euler);
+            //printf("  ** done refinements **\n");
+
+            //printf("maps after refinement\n");
+            //show_maps(mesh,coarse_map,fine_map);
+         }
+
+
+         cur_ref_set = new_ref_set;
+
+         // Perform any new derefinements
+         if (coarsen_set.size()) {
+
+            //printf("  ** start derefinements **\n");
+
+            Array<double> mock_error(mesh.GetNE());
+            mock_error = 1.0;
+
+            // We only mark one of the fine elements, but that's fine
+            // because we can set the threshold low enough to always
+            // derefine.
+            for (it = coarsen_set.begin(); it != coarsen_set.end(); ++it) {
+               const int& i = *it;
+               //printf("coarsening ref idx %d\n",i);
+               assert(coarse_map[i] == -1);
+               int i1 = fine_map[i];
+               //printf("which is fine element %d (and +1)\n",i1);
+               assert(i1 >= 0);
+               mock_error[i1] = 0.0;
+               //printf("setting mock error to 0.0 in %d\n",i1);
+            }
+            mesh.DerefineByError(mock_error, 2.0);
+            amr_update(fes, dfes, vfes, u_block, rho, rho_u, rho_e, sol);
+            A.Update();
+            Aflux.Update();
+            Aflux.Assemble();
+            euler.Update();
+            ode_solver->Init(euler);
+
+            //printf("updating coarse_map and fine_map after derefinement\n");
+
+            const CoarseFineTransformations& cft = mesh.ncmesh->GetDerefinementTransforms();
+
+            Table c2f;
+            cft.GetCoarseToFineMap(mesh, c2f);
+            //c2f.Print();
+
+            Array<int> row;
+            map<int,int> old2new;
+            for (int i = 0; i < c2f.Size(); ++i) {
+               c2f.GetRow(i,row);
+               for (int j = 0; j < row.Size(); j++) {
+                  old2new[row[j]] = i;
+               }
+            }
+
+            // map<int,int>::iterator it;
+            // printf("old2new\n");
+            // for (it = old2new.begin(); it != old2new.end(); ++it) {
+            //    printf("%d -> %d\n",it->first,it->second);
+            // }
+
+            for (size_t i = 0; i < coarse_map.size(); i++) {
+               //printf("ref i = %lu\n",i);
+               int jfine = fine_map[i];
+               int jcoarse = coarse_map[i];
+               //printf("mesh j old grid (fine) = %d\n",jfine);
+               if (jfine > -1) {
+                  int newj = old2new[jfine];
+                  //printf("mesh j new grid (coarse) = %d\n",newj);
+                  c2f.GetRow(newj,row);
+                  if (row.Size() > 1) {
+                     //printf("  was fine, now coarse\n");
+                     // was fine, now coarse
+                     coarse_map[i] = newj;
+                     fine_map[i] = -1;
+                  }
+                  else {
+                     //printf("  was fine, still fine\n");
+                     // was fine, still fine
+                     fine_map[i] = newj;
+                  }
+               }
+               //printf("mesh j old grid coarse = %d\n",jcoarse);
+               if (jcoarse > -1) {
+                  int newj = old2new[jcoarse];
+                  coarse_map[i] = newj;
+               }
+            }
+
+            // printf("new coarse_map is:\n");
+            // for (size_t i = 0; i < coarse_map.size(); i++) {
+            //    printf("%lu -> %d\n",i,coarse_map[i]);
+            // }
+            // printf("new fine_map is:\n");
+            // for (size_t i = 0; i < fine_map.size(); i++) {
+            //    printf("%lu -> %d\n",i,fine_map[i]);
+            // }
+
+            //printf("  ** done derefinements **\n");
+
+            //printf("maps after derefinement\n");
+            //show_maps(mesh,coarse_map,fine_map);
+
+         }
+
+      }
+
+      // Output the current solution.
+      if (output_cycle_soln) {
+
+         ostringstream fn;
+         fn << "soln-" << ti << ".mesh";
+         ofstream mesh_ofs(fn.str());
+         mesh_ofs.precision(precision);
+         mesh_ofs << mesh;
+
+         for (int k = 0; k < num_equation; k++)
+         {
+            GridFunction uk(&fes, u_block.GetBlock(k));
+            ostringstream sol_name;
+            sol_name << "soln-" << ti << "-" << k << ".gf";
+            ofstream sol_ofs(sol_name.str().c_str());
+            sol_ofs.precision(precision);
+            sol_ofs << uk;
+         }
+      }
+
+      // Output the current error on the reference mesh
+      if (output_cycle_errors) {
+         vector<double> ref_errors(coarse_map.size());
+         compute_reference_errors(ti, fec, mesh, rho,
+                                  coarse_map, fine_map, ref_errors);
+         ostringstream fn;
+         fn << "err-" << ti << ".dat";
+         ofstream err_ofs(fn.str());
+         vector<double>::iterator it;
+         for (it = ref_errors.begin(); it != ref_errors.end(); ++it) {
+            int i = std::distance(ref_errors.begin(), it);
+            double err = *it;
+            // "mark" refinements with negation. We'll use this to vis
+            // the refined regions on the error plots.
+            if (coarse_map[i] < 0) err *= -1;
+            err_ofs << i << " " << err << endl;
+         }
+      }
+
+
+      if (cfl > 0) {
+         dt = estimate_dt(cfl, mesh, sol, A, order);
+      }
+
       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);
@@ -275,9 +774,10 @@ int main(int argc, char *argv[])
          cout << "time step: " << ti << ", time: " << t << endl;
          if (visualization)
          {
-            sout << "solution\n" << mesh << mom << flush;
+            sout << "solution\n" << mesh << rho_u << flush;
          }
       }
+
    }
 
    tic_toc.Stop();
@@ -285,21 +785,57 @@ int main(int argc, char *argv[])
 
    // 9. Save the final solution. This output can be viewed later using GLVis:
    //    "glvis -m vortex.mesh -g vortex-1-final.gf".
-   for (int k = 0; k < num_equation; k++)
    {
-      GridFunction uk(&fes, u_block.GetBlock(k));
-      ostringstream sol_name;
-      sol_name << "vortex-" << k << "-final.gf";
-      ofstream sol_ofs(sol_name.str().c_str());
-      sol_ofs.precision(precision);
-      sol_ofs << uk;
+      ofstream mesh_ofs("vortex-final.mesh");
+      mesh_ofs.precision(precision);
+      mesh_ofs << mesh;
+
+
+      for (int k = 0; k < num_equation; k++)
+      {
+         GridFunction uk(&fes, u_block.GetBlock(k));
+         ostringstream sol_name;
+         sol_name << "vortex-" << k << "-final.gf";
+         ofstream sol_ofs(sol_name.str().c_str());
+         sol_ofs.precision(precision);
+         sol_ofs << uk;
+      }
    }
 
    // 10. Compute the L2 solution error summed for all components.
-   if (t_final == 2.0)
+   if (rseq.Size())
    {
-      const double error = sol.ComputeLpError(2, u0);
-      cout << "Solution error: " << error << endl;
+      // const CoarseFineTransformations& cft = mesh.GetRefinementTransforms();
+      // Table c2f;
+      // cft.GetCoarseToFineMap(mesh, c2f);
+
+      // Refine to the reference mesh by refining every coarse el that
+      // wasn't refined.
+      Array<int> refs;
+      for (int i = 0; i < mesh.GetNE(); i++) {
+         if (mesh.ncmesh->GetElementDepth(i) == 0) {
+            refs.Append(i);
+         }
+      }
+      mesh.GeneralRefinement(refs);
+
+      amr_update(fes, dfes, vfes, u_block, rho, rho_u, rho_e, sol);
+
+      Mesh mesh_ref("reference.mesh");
+      ifstream rho_ifs("reference-rho.gf");
+      ifstream rho_u_ifs("reference-rho-u.gf");
+      ifstream rho_e_ifs("reference-rho-e.gf");
+      GridFunction rho_ref(&mesh_ref, rho_ifs);
+      GridFunction rho_u_ref(&mesh_ref, rho_u_ifs);
+      GridFunction rho_e_ref(&mesh_ref, rho_e_ifs);
+      rho_ref.Print();
+      rho.Print();
+
+      GridFunctionCoefficient rho_ref_coeff(&rho_ref);
+      const double err = rho.ComputeL2Error(rho_ref_coeff);
+      printf("final L2 error %e\n",err);
+      // const double error = sol.ComputeLpError(2, u0);
+      // cout << "Solution error: " << error << endl;
    }
 
    // Free the used memory.

examples/ex18.hpp

@@ -11,7 +11,8 @@ extern int problem;
 // Maximum characteristic speed (updated by integrators)
 extern double max_char_speed;
 
-extern const int num_equation;
+extern int num_equation;
+
 extern const double specific_heat_ratio;
 extern const double gas_constant;
 
@@ -24,7 +25,8 @@ private:
 
    FiniteElementSpace &vfes;
    Operator &A;
-   SparseMatrix &Aflux;
+   MixedBilinearForm& Aflux;
+   //SparseMatrix &Aflux;
    DenseTensor Me_inv;
 
    mutable Vector state;
@@ -36,10 +38,13 @@ private:
 
 public:
    FE_Evolution(FiniteElementSpace &vfes_,
-                Operator &A_, SparseMatrix &Aflux_);
+                Operator &A_,
+                MixedBilinearForm &Aflux_);
 
    virtual void Mult(const Vector &x, Vector &y) const;
 
+   void Update();
+
    virtual ~FE_Evolution() { }
 };
 
@@ -100,7 +105,7 @@ public:
 
 // Implementation of class FE_Evolution
 FE_Evolution::FE_Evolution(FiniteElementSpace &vfes_,
-                           Operator &A_, SparseMatrix &Aflux_)
+                           Operator &A_, MixedBilinearForm &Aflux_)
    : TimeDependentOperator(A_.Height()),
      dim(vfes_.GetFE(0)->GetDim()),
      vfes(vfes_),
@@ -123,6 +128,32 @@ FE_Evolution::FE_Evolution(FiniteElementSpace &vfes_,
       inv.Factor();
       inv.GetInverseMatrix(Me_inv(i));
    }
+   //Aflux.Assemble();
+}
+
+void FE_Evolution::Update()
+{
+   height = width = A.Height();
+
+   Me_inv.SetSize(vfes.GetFE(0)->GetDof(),
+                  vfes.GetFE(0)->GetDof(),
+                  vfes.GetNE());
+
+   Aflux.Assemble();
+
+   // Standard local assembly and inversion for energy mass matrices.
+   const int dof = vfes.GetFE(0)->GetDof();
+   DenseMatrix Me(dof);
+   DenseMatrixInverse inv(&Me);
+   MassIntegrator mi;
+   for (int i = 0; i < vfes.GetNE(); i++)
+   {
+      mi.AssembleElementMatrix(*vfes.GetFE(i), *vfes.GetElementTransformation(i), Me);
+      inv.Factor();
+      inv.GetInverseMatrix(Me_inv(i));
+   }
+
+   z.SetSize(A.Height());
 }
 
 void FE_Evolution::Mult(const Vector &x, Vector &y) const
@@ -385,6 +416,8 @@ void FaceIntegrator::AssembleFaceVector(const FiniteElement &el1,
                                         FaceElementTransformations &Tr,
                                         const Vector &elfun, Vector &elvect)
 {
+  int dim = el1.GetDim();
+
    // Compute the term <F.n(u),[w]> on the interior faces.
    const int dof1 = el1.GetDof();
    const int dof2 = el2.GetDof();
@@ -424,16 +457,25 @@ void FaceIntegrator::AssembleFaceVector(const FiniteElement &el1,
 
       Tr.SetAllIntPoints(&ip); // set face and element int. points
 
+      const IntegrationPoint &eip1 = Tr.GetElement1IntPoint();
+      const IntegrationPoint &eip2 = Tr.GetElement2IntPoint();
+
       // Calculate basis functions on both elements at the face
-      el1.CalcShape(Tr.GetElement1IntPoint(), shape1);
-      el2.CalcShape(Tr.GetElement2IntPoint(), shape2);
+      el1.CalcShape(eip1, shape1);
+      el2.CalcShape(eip2, 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);
+      if (1 == dim) {
+	nor(0) = 2*eip1.x - 1.0;
+      }
+      else {
+	CalcOrtho(Tr.Jacobian(), nor);
+      }
+
       const double mcs = rsolver.Eval(funval1, funval2, nor, fluxN);
 
       // Update max char speed
@@ -504,62 +546,147 @@ bool StateIsPhysical(const Vector &state, const int dim)
 // Initial condition
 void InitialCondition(const Vector &x, Vector &y)
 {
-   MFEM_ASSERT(x.Size() == 2, "");
+  int dim = x.Size();
 
-   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");
-   }
+  if (dim == 1) {
 
-   const double xc = 0.0, yc = 0.0;
+     double rho;
+     double u;
+     double e;
 
-   // Nice units
-   const double vel_inf = 1.;
-   const double den_inf = 1.;
+     if (problem == 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 xc = x[0]-0.5;
+        rho = 1.0;
+        u = 0.0;
+        e = 1.0+0.1*exp(-xc*xc);
 
-   double r2rad = 0.0;
-   r2rad += (x(0) - xc) * (x(0) - xc);
-   r2rad += (x(1) - yc) * (x(1) - yc);
-   r2rad /= (radius * radius);
+     }
+     else if (problem == 2) {
 
-   const double shrinv1 = 1.0 / (specific_heat_ratio - 1.);
+        u = 0.0;
+        double gamma = 1.4;
+        double gm1 = gamma -1.0;
 
-   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;
+        double rhoL = 1.0;
+        double pL = 1.0;
 
-   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);
+        double rhoR = 0.125;
+        double pR = 0.1;
 
-   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;
+        if (x[0] < 0.5) {
+           rho = rhoL;
+           e = 1./gm1*pL/rhoL;
+        }
+        else {
+           rho = rhoR;
+           e = 1./gm1*pR/rhoR;
+        }
+     }
+     else if (problem == 3) {
+
+        u = 0.0;
+        double gamma = 1.4;
+        double gm1 = gamma -1.0;
+
+        double rhoL = 1.0;
+        double pL = 1.0;
+
+        double rhoR = 0.125;
+        double pR = 0.1;
+
+        rho = rhoR;
+        e = 1./gm1*pR/rhoR;
+        if (x[0] > 0.3 && x[0] < 0.4) {
+           rho = rhoL;
+           e = 1./gm1*pL/rhoL;
+        }
+        if (x[0] > 0.5 && x[0] < 1.0) {
+           double xc = 0.75;
+           double x0 = x[0]-xc;
+           double a = 10.0;
+           double x1 = a*x0;
+           rho = rhoR +0.5*rhoL*exp(-x1*x1);
+           e = 1./gm1*pR/rhoL;
+        }
+     }
+     else {
+        printf("unknown problem!\n");
+        exit(1);
+     }
+
+    // if (x(0) < 0.5) {
+    //   rho = 1.0;
+    //   e = 1.05*(;
+    // }
+
+    int neq = 0;
+    y(neq++) = rho;
+    y(neq++) = rho * u;
+    y(neq++) = rho * e;
+  }
+  else {
+
+    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;
+
+    int neq = 0;
+    y(neq++) = den;
+    y(neq++) = den * velX;
+    y(neq++) = den * velY;
+    y(neq++) = den * energy;
+
+  } // 2d
 
-   y(0) = den;
-   y(1) = den * velX;
-   y(2) = den * velY;
-   y(3) = den * energy;
 }

examples/generate-glvis-script

@@ -0,0 +1,48 @@
+#!/bin/bash
+
+if [ $# -ne 2 ]; then
+    echo "Usage: generate-glvis-script <ncyc> <output-file-prefix>"
+fi
+
+ncyc=$1
+prefix=$2
+comp=0
+
+ncycm1=$(( $ncyc - 1 ))
+
+cat <<EOF > glvis-script
+window 0 0 400 400
+
+# Initial solution
+solution soln-$comp.mesh soln-0-$comp.gf
+
+# Setup the GLVis scene. Executed after pressing the space bar.
+{
+   perspective off
+   viewcenter 0.2 0
+   valuerange 0.0 1.0
+   keys ca
+   autoscale off
+}
+
+# Take multiple screenshots. Executed after pressing space bar.
+{
+EOF
+
+for i in `seq 0 $ncycm1`; do
+    if [ ! -f soln-$i.mesh ]; then
+	echo "The file soln-$i.mesh doesn't exist. Stopping."
+	exit 1
+    fi
+    if [ ! -f soln-$i-$comp.gf ]; then
+	echo "The file soln-$i-$comp.gf doesn't exist. Stopping."
+	exit 1
+    fi
+    ii=`printf %04d $i`
+    echo "solution soln-$i.mesh soln-$i-$comp.gf" >> glvis-script
+    #echo "valuerange 0.0 1.0" >> glvis-script
+    echo "screenshot $prefix-soln-$ii.png" >> glvis-script
+done
+echo "}" >> glvis-script
+
+echo "% glvis -run glvis-script"

examples/generate-screenshot-cmds

@@ -0,0 +1,11 @@
+#!/bin/bash
+
+rm -f tmp
+for i in {0..19}; do
+    ii=`printf %04d $i`
+    echo "solution soln-$i.mesh soln-$i-0.gf" >> tmp
+    echo "valuerange 0.0 1.0" >> tmp
+    echo "screenshot greedy-soln-$ii.png" >> tmp
+done
+
+echo "output in ./tmp"

examples/make-err-plots.py

@@ -0,0 +1,32 @@
+
+# Take err-N.dat files and create error plots for each cycle.
+# Negative values are interpreted as refined regions and plotted in
+# red.
+
+import math
+import matplotlib.pyplot as plt
+
+for i in range(20):
+    print(f"{i}")
+    with open(f"err-{i}.dat") as f:
+        lines = f.readlines()
+        x = [float(line.split()[0]) for line in lines]
+        y = [float(line.split()[1]) for line in lines]
+
+        # errors in refined zones are marked as negative by convention
+        # from ex18. in order to detect "negative zero", we use the
+        # copysign method, which distinguishes negative zero from
+        # zero.
+        r = [(x,-y) for x,y in zip(x,y) if math.copysign(1, y) == -1.0]
+        xr,yr = zip(*r)
+        print(xr)
+        print(yr)
+        plt.scatter(x,y)
+        plt.scatter(xr,yr,c='r')
+        plt.title(f"Timestep {i}")
+        plt.xlabel("reference el")
+        plt.ylabel("L2 error")
+        plt.ylim(0,0.05)
+        #plt.show()
+        plt.savefig(f"opt-err-{i}.png",dpi=60)
+        plt.clf()

examples/make-glvis-figs

@@ -0,0 +1,19 @@
+window 0 0 400 400
+
+# Initial solution
+solution soln-0.mesh soln-0-0.gf
+
+# Setup the GLVis scene. Executed after pressing the space bar.
+{
+   perspective off
+#   view 0 0
+   viewcenter 0.2 0
+   valuerange 0.0 1.0
+#   zoom 1.5
+   keys ca
+}
+
+# Take multiple screenshots. Executed after pressing the space bar.
+{
+#insert the screenshot commands here
+}

examples/make-periodic-mesh

@@ -0,0 +1,63 @@
+#!/bin/bash
+
+if [ $# -ne 1 ]; then
+    echo "Usage: ./make-periodic-mesh <NE>"
+fi
+
+NE=$1
+NEm1=$(( $NE - 1 ))
+NEm2=$(( $NE - 2 ))
+
+filename="periodic-segment-$NE.mesh"
+
+cat <<ENDHERE > $filename
+MFEM mesh v1.0
+
+#
+# MFEM Geometry Types (see mesh/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+#
+
+dimension
+1
+
+# format: <attribute> <geometry type> <vertex 0> <vertex 1> ...
+elements
+ENDHERE
+
+echo "$NE" >> $filename
+for i in `seq 0 $NEm2`; do
+    echo "1 1 $i $((i+1))" >> $filename
+done
+echo "1 1 $((i+1)) 0" >> $filename
+
+cat <<EOF >> $filename
+
+boundary
+0
+
+vertices
+$NE
+
+nodes
+FiniteElementSpace
+FiniteElementCollection: L2_T1_1D_P1
+VDim: 1
+Ordering: 1
+ 
+EOF
+
+dx=`echo "1.0 / $NE" | bc -l`
+for i in `seq 0 $NEm1`; do
+    p1=`echo "$i * $dx" | bc -l`
+    p2=`echo "$((i+1)) * $dx" | bc -l`
+    echo "$p1 $p2" >> $filename
+done
+
+

examples/periodic-segment-16.mesh

@@ -0,0 +1,65 @@
+MFEM mesh v1.0
+
+#
+# MFEM Geometry Types (see mesh/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+#
+
+dimension
+1
+
+# format: <attribute> <geometry type> <vertex 0> <vertex 1> ...
+elements
+16
+1 1 0 1
+1 1 1 2
+1 1 2 3
+1 1 3 4
+1 1 4 5
+1 1 5 6
+1 1 6 7
+1 1 7 8
+1 1 8 9
+1 1 9 10
+1 1 10 11
+1 1 11 12
+1 1 12 13
+1 1 13 14
+1 1 14 15
+1 1 15 0
+
+
+boundary
+0
+
+vertices
+16
+
+nodes
+FiniteElementSpace
+FiniteElementCollection: L2_T1_1D_P1
+VDim: 1
+Ordering: 1
+
+0.0000 0.0625
+0.0625 0.1250
+0.1250 0.1875
+0.1875 0.2500
+0.2500 0.3125
+0.3125 0.3750
+0.3750 0.4375
+0.4375 0.5000
+0.5000 0.5625
+0.5625 0.6250
+0.6250 0.6875
+0.6875 0.7500
+0.7500 0.8125
+0.8125 0.8750
+0.8750 0.9375
+0.9375 1.0000

examples/periodic-segment-8.mesh

@@ -0,0 +1,50 @@
+MFEM mesh v1.0
+
+#
+# MFEM Geometry Types (see mesh/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+#
+
+dimension
+1
+
+# format: <attribute> <geometry type> <vertex 0> <vertex 1> ...
+elements
+8
+1 1 0 1
+1 1 1 2
+1 1 2 3
+1 1 3 4
+1 1 4 5
+1 1 5 6
+1 1 6 7
+1 1 7 0
+
+boundary
+0
+
+vertices
+8
+
+nodes
+FiniteElementSpace
+FiniteElementCollection: L2_T1_1D_P1
+VDim: 1
+Ordering: 1
+
+0.000 0.125
+0.125 0.250
+0.250 0.375
+0.375 0.500
+0.500 0.625
+0.625 0.750
+0.750 0.875
+0.875 1.000
+
+

examples/run-search

@@ -0,0 +1,43 @@
+#!/bin/bash
+
+# sod 5 regrids
+options="-dt 0.01 --cfl-number -1 " # fixed timestep
+options+="--regrid-period 4 "       # regrid period
+options+="--problem 2 "             # Sod problem
+options+="--order 0 "               # 0th order DG
+options+="--refine 0 "              # no initial uniform refinement
+options+="--t-final 0.20 "
+
+# create reference solution on 16 element mesh
+./ex18 $options -m periodic-segment-20.mesh
+./save-reference
+
+ne=10
+nem1=$((ne-1))
+
+echo "" > err-results.dat
+for i1 in $(seq 0 $nem1); do
+    for i2 in $(seq 0 $nem1); do
+	for i3 in $(seq 0 $nem1); do
+	    for i4 in $(seq 0 $nem1); do
+		for i5 in $(seq 0 $nem1); do
+		    dir="seq-$i1$i2$i3$i4$i5"
+		    mkdir $dir
+		    cd $dir
+		    ln -s ../reference.mesh .
+		    ln -s ../reference-rho.gf .
+		    ln -s ../reference-rho-u.gf .
+		    ln -s ../reference-rho-e.gf .
+		    echo "running $i1$i2$i3$i4$i5"
+		    seq="$i1 $i2 $i3 $i4 $i5"
+		    cmd="../ex18 -rseq \"$seq\" $options -m ../periodic-segment-10.mesh | grep err | awk \"{print \$2}\""
+		    echo "$cmd" > runcode
+		    err=$( source runcode )
+		    cd ..
+		    echo "$i1 $i2 $i3 $i4 $i5 $err" >> err-results.dat
+		done
+	    done
+	done
+    done
+done
+

examples/save-reference

@@ -0,0 +1,16 @@
+#!/bin/bash
+
+# The final time is the reference solution
+mv vortex.mesh reference.mesh
+mv vortex-0-final.gf reference-rho.gf
+mv vortex-1-final.gf reference-rho-u.gf
+mv vortex-2-final.gf reference-rho-e.gf
+
+# Save the timesteps as well for errors
+for i in {0..19}
+do
+    mv soln-$i.mesh reference-$i.mesh
+    mv soln-$i-0.gf reference-$i-0.gf
+    mv soln-$i-1.gf reference-$i-1.gf
+    mv soln-$i-2.gf reference-$i-2.gf
+done

mesh/mesh.cpp

@@ -3495,8 +3495,6 @@ void Mesh::SetMeshGen()
             mesh_geoms |= (1 << Geometry::TETRAHEDRON);
          case Element::TRIANGLE:
             mesh_geoms |= (1 << Geometry::TRIANGLE);
-         case Element::SEGMENT:
-            mesh_geoms |= (1 << Geometry::SEGMENT);
          case Element::POINT:
             mesh_geoms |= (1 << Geometry::POINT);
             meshgen |= 1;
@@ -3508,6 +3506,8 @@ void Mesh::SetMeshGen()
             mesh_geoms |= (1 << Geometry::SQUARE);
             mesh_geoms |= (1 << Geometry::SEGMENT);
             mesh_geoms |= (1 << Geometry::POINT);
+         case Element::SEGMENT:
+            mesh_geoms |= (1 << Geometry::SEGMENT);
             meshgen |= 2;
             break;
 

mesh/ncmesh.cpp

@@ -4253,6 +4253,19 @@ void NCMesh::GetPointMatrix(Geometry::Type geom, const char* ref_path,
             pm = PointMatrix(mid12, mid20, mid01);
          }
       }
+      else if (geom == Geometry::SEGMENT)
+      {
+         Point mid01(pm(0), pm(1));
+
+         if (child == 0)
+         {
+            pm = PointMatrix(pm(0), mid01);
+         }
+         else if (child == 1)
+         {
+            pm = PointMatrix(mid01, pm(1));
+         }
+      }
    }
 
    // write the points to the matrix