merge with zz_test branch

Simplify Mesh::GetFaceElements

ZZ_test/ZZ_test.cpp

@@ -0,0 +1,409 @@
+//                         Test of ZZ error estimator
+//
+// Compile with: make ZZ_test
+//
+
+#include "mfem.hpp"
+#include <fstream>
+
+// #include "exact.hpp"
+
+using namespace std;
+using namespace mfem;
+
+double lshape_exsol(const Vector &p);
+void   lshape_exgrad(const Vector &p, Vector &grad);
+double lshape_laplace(const Vector &p);
+
+double sinsin_exsol(const Vector &p);
+void   sinsin_exgrad(const Vector &p, Vector &grad);
+double sinsin_laplace(const Vector &p);
+
+double poly_exsol(const Vector &p);
+void   poly_exgrad(const Vector &p, Vector &grad);
+double poly_laplace(const Vector &p);
+
+int dim;
+const char* keys = "Rjlmc*******";
+
+int main(int argc, char *argv[])
+{
+   // Parse command-line options.
+   int problem = 0;
+   int order = 1;
+   double ref_threshold = 0.8;
+//    const char *elemerr_file = "elemerr.txt";
+   int nc_limit = 1;
+   const char *device_config = "cpu";
+   bool visualization = false;
+   int which_estimator = 0;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&problem, "-p", "--problem",
+                  "Problem type: 0 = canonical L-shaped solution, 1 = sinusoid.");
+   args.AddOption(&order, "-o", "--order",
+                  "Initial mesh finite element order (polynomial degree).");
+   args.AddOption(&ref_threshold, "-rt", "--ref-threshold",
+                  "Refine elements with error larger than threshold * max_error.");
+   args.AddOption(&nc_limit, "-nc", "--nc-limit",
+                  "Set maximum difference of refinement levels of adjacent elements.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&which_estimator, "-est", "--estimator",
+                  "Which estimator to use: "
+                  "0 = ZZ, 1 = Kelly. Defaults to ZZ.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   // Enable hardware devices such as GPUs, and programming models such as
+   // CUDA, OCCA, RAJA and OpenMP based on command line options.
+   Device device(device_config);
+   device.Print();
+
+   const char *mesh_file;
+   // if (problem == 0)
+   // {
+      mesh_file = "l-shape-benchmark.mesh";
+   // }
+
+   // 2. Read the (serial) mesh from the given mesh file.
+   Mesh mesh(mesh_file, 1, 1);
+   dim = mesh.Dimension();
+   mesh.EnsureNCMesh();
+   mesh.UniformRefinement(); // ZZ doesn't work properly on the initial L-shaped mesh
+
+   // Define a finite element space on the mesh.
+   H1_FECollection fec(order, dim);
+   L2_FECollection l2fec(order, dim);
+   FiniteElementSpace fespace(&mesh, &fec);
+
+   // Define the solution vector x as a finite element grid function
+   // corresponding to fespace.
+   GridFunction x(&fespace);
+
+   // Define exact solutions
+   FunctionCoefficient *exsol=nullptr;
+   VectorFunctionCoefficient *exgrad=nullptr;
+   FunctionCoefficient *rhs=nullptr;
+   ConstantCoefficient one(1.0);
+
+   switch (problem)
+   {
+      case 1:
+      {
+         exsol = new FunctionCoefficient(sinsin_exsol);
+         exgrad = new VectorFunctionCoefficient(dim, sinsin_exgrad);
+         rhs = new FunctionCoefficient(sinsin_laplace);
+         break;
+      }
+      case 2:
+      {
+         exsol = new FunctionCoefficient(poly_exsol);
+         exgrad = new VectorFunctionCoefficient(dim, poly_exgrad);
+         rhs = new FunctionCoefficient(poly_laplace);
+         break;
+      }
+      default:
+      case 0:
+      {
+         exsol = new FunctionCoefficient(lshape_exsol);
+         exgrad = new VectorFunctionCoefficient(dim, lshape_exgrad);
+         rhs = new FunctionCoefficient(lshape_laplace);
+         break;
+      }
+   }
+
+   // Set up the linear form b(.) and the bilinear form a(.,.).
+   LinearForm b(&fespace);
+   BilinearForm a(&fespace);
+
+   b.AddDomainIntegrator(new DomainLFIntegrator(*rhs));
+   a.AddDomainIntegrator(new DiffusionIntegrator(one));
+
+   // All boundary attributes will be used for essential (Dirichlet) BC.
+   MFEM_VERIFY(mesh.bdr_attributes.Size() > 0,
+               "Boundary attributes required in the mesh.");   
+
+   // Connect to GLVis.
+   socketstream sol_sock, ord_sock, dbg_sock[3], err_sock;
+
+   ostringstream file_name;
+   file_name << "conv_order" << order << ".csv";
+   ofstream conv(file_name.str().c_str());
+//    std::ofstream elemerr(elemerr_file);
+
+   if (visualization)
+   {
+      cout << "\n Press enter to advance... " << endl;
+   }
+
+   cout << setw(4) << "\nRef." << setw(12) << "DOFs" << setw(21) << "H^1_0 error" << setw(21) << "error estimate" << setw(18) << "H^1_0 rate" << setw(18) << "estimator rate" <<  endl;
+   conv << "DOFs " << ", " << "H^1_0 error" << ", " << "error estimate" << endl;
+
+   double old_num_dofs = 0.0;
+   double old_H10_error = 0.0;
+   double old_ZZ_error = 0.0;
+   const int max_dofs = 20000;
+   for (int it = 0; ; it++)
+   {
+      int num_dofs = fespace.GetTrueVSize();
+
+      // Set Dirichlet boundary values in the GridFunction x.
+      // Determine the list of Dirichlet true DOFs in the linear system.
+      Array<int> ess_bdr(mesh.bdr_attributes.Max());
+      ess_bdr = 1;
+      x = 0.0;
+      x.ProjectBdrCoefficient(*exsol, ess_bdr);
+      Array<int> ess_tdof_list;
+      fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   
+      // Solve for the current mesh:
+      b.Assemble();
+      a.Assemble();
+      a.Finalize();
+      OperatorPtr A;
+      Vector B, X;
+   
+      const int copy_interior = 1;
+      a.FormLinearSystem(ess_tdof_list, x, b, A, X, B, copy_interior);
+      GSSmoother M((SparseMatrix&)(*A));
+      PCG(*A, M, B, X, 0, 2000, 1e-30, 0.0);
+
+      a.RecoverFEMSolution(X, b, x);
+
+      // Calculate the total error in the H^1_0 norm.
+      double H10_error = x.ComputeGradError(exgrad);
+      DiffusionIntegrator di;
+
+      ErrorEstimator* estimator{nullptr};
+      switch (which_estimator)
+      {
+         case 1:
+         {
+            auto flux_fes = new FiniteElementSpace(&mesh, &fec, dim);
+            estimator = new NewZienkiewiczZhuEstimator(di, x, flux_fes);
+            // int flux_order = 4;
+            // estimator = new NewZienkiewiczZhuEstimator(di, x, flux_order);
+            break;
+         }
+
+         case 2:
+         {
+            auto flux_fes = new FiniteElementSpace(&mesh, &l2fec, dim);
+            estimator = new KellyErrorEstimator(di, x, flux_fes);
+            break;
+         }
+
+         default:
+            std::cout << "Unknown estimator. Falling back to ZZ." << std::endl;
+         case 0:
+         {
+            auto flux_fes = new FiniteElementSpace(&mesh, &fec, dim);
+            estimator = new ZienkiewiczZhuEstimator(di, x, flux_fes);
+            break;
+         }
+      }
+      StopWatch chrono;
+      chrono.Clear();
+      chrono.Start();
+      const Vector &zzerr = estimator->GetLocalErrors();
+      chrono.Stop();
+      mfem::out << "get local errors time: " << chrono.RealTime() << endl;
+      // double ZZ_error = zzerr.Norml2();
+      double ZZ_error = estimator->GetTotalError();
+
+      // estimate convergence rate
+      double H10_rate = 0.0;
+      double ZZ_rate = 0.0;
+      if (old_H10_error > 0.0)
+      {
+          H10_rate = log(H10_error/old_H10_error) / log(old_num_dofs/num_dofs);
+          ZZ_rate  = log(ZZ_error/old_ZZ_error)   / log(old_num_dofs/num_dofs);
+      }
+
+      cout << setw(4) << it << setw(12) << num_dofs << setw(21) << H10_error << setw(21) << ZZ_error << setw(18) << H10_rate << setw(18) << ZZ_rate << endl;
+
+      // Send solution by socket to the GLVis server.
+      if (visualization)
+      {
+         cin.get();
+         const char vishost[] = "localhost";
+         const int  visport   = 19916;
+         socketstream sol_sock(vishost, visport);
+         sol_sock.precision(8);
+         sol_sock << "solution\n" << mesh << x;
+         sol_sock << "keys ARjlm\n";
+      }
+
+      if (num_dofs > max_dofs)
+      {
+         cout << "\n Reached the maximum number of dofs. Stop." << endl;
+         break;
+      }
+
+      // Save dofs and error for convergence plot
+      conv << num_dofs << ", " << H10_error << ", " << ZZ_error << endl;
+
+    //   for (int i = 0; i < mesh.GetNE(); i++)
+    //   {
+    //      elemerr << sqrt(elemError[i]) << ' ';
+    //   }
+    //   elemerr << endl;
+
+
+      Array<Refinement> refinements;
+      double err_max = zzerr.Max();
+      for (int i = 0; i < mesh.GetNE(); i++)
+      {
+         if (zzerr[i] > ref_threshold * err_max)
+         {
+             refinements.Append(Refinement(i, 7));
+         }
+      }
+      mesh.GeneralRefinement(refinements, -1, nc_limit);
+
+      old_num_dofs = double(num_dofs);
+      old_H10_error = H10_error;
+      old_ZZ_error = ZZ_error;
+
+      // Update the space, interpolate the solution.
+      fespace.Update();
+      a.Update();
+      b.Update();
+      x.Update();
+
+      // Free the used memory.
+      delete estimator;
+   }
+
+   // Free the used memory.
+   delete exsol;
+   delete exgrad;
+   delete rhs;
+   return 0;
+}
+
+
+// L-shape domain problem exact solution (2D)
+
+double lshape_exsol(const Vector &p)
+{
+   double x = p(0), y = p(1);
+   double r = sqrt(x*x + y*y);
+   double a = atan2(y, x);
+   if (a < 0) { a += 2*M_PI; }
+   return pow(r, 2.0/3.0) * sin(2.0*a/3.0);
+}
+
+void lshape_exgrad(const Vector &p, Vector &grad)
+{
+   double x = p(0), y = p(1);
+   double a = atan2(y, x);
+   if (a < 0) { a += 2*M_PI; }
+   double theta23 = 2.0/3.0*a;
+   double r23 = pow(x*x + y*y, 2.0/3.0);
+   grad(0) = 2.0/3.0*x*sin(theta23)/(r23) - 2.0/3.0*y*cos(theta23)/(r23);
+   grad(1) = 2.0/3.0*y*sin(theta23)/(r23) + 2.0/3.0*x*cos(theta23)/(r23);
+}
+
+double lshape_laplace(const Vector &p)
+{
+   return 0;
+}
+
+double sinsin_exsol(const Vector &p)
+{
+   double x = p(0), y = p(1);
+   return sin(M_PI * x) * sin(M_PI * y);
+}
+
+void sinsin_exgrad(const Vector &p, Vector &grad)
+{
+   double x = p(0), y = p(1);
+   grad(0) = M_PI * cos(M_PI * x) * sin(M_PI * y);
+   grad(1) = M_PI * sin(M_PI * x) * cos(M_PI * y);
+}
+
+double sinsin_laplace(const Vector &p)
+{
+   double x = p(0), y = p(1);
+   return 2 * M_PI * M_PI * sin(M_PI * x) * sin(M_PI * y);
+}
+
+// double poly_exsol(const Vector &p)
+// {
+//    double x = p(0), y = p(1);
+//    return y*y;
+// }
+
+// void poly_exgrad(const Vector &p, Vector &grad)
+// {
+//    double x = p(0), y = p(1);
+//    grad(0) = 0.0;
+//    grad(1) = 2.0*y;
+// }
+
+// double poly_laplace(const Vector &p)
+// {
+//    return -2.0;
+// }
+
+// double poly_exsol(const Vector &p)
+// {
+//    double x = p(0), y = p(1);
+//    return x*x;
+// }
+
+// void poly_exgrad(const Vector &p, Vector &grad)
+// {
+//    double x = p(0), y = p(1);
+//    grad(0) = 2.0*x;
+//    grad(1) = 0.0;
+// }
+
+// double poly_laplace(const Vector &p)
+// {
+//    return -2.0;
+// }
+
+double poly_exsol(const Vector &p)
+{
+   double x = p(0), y = p(1);
+   return x*y;
+}
+
+void poly_exgrad(const Vector &p, Vector &grad)
+{
+   double x = p(0), y = p(1);
+   grad(0) = y;
+   grad(1) = x;
+}
+
+double poly_laplace(const Vector &p)
+{
+   return 0.0;
+}
+
+// double poly_exsol(const Vector &p)
+// {
+//    double x = p(0), y = p(1);
+//    return x + 2.0*y;
+// }
+
+// void poly_exgrad(const Vector &p, Vector &grad)
+// {
+//    double x = p(0), y = p(1);
+//    grad(0) = 1.0;
+//    grad(1) = 2.0;
+// }
+
+// double poly_laplace(const Vector &p)
+// {
+//    return 0.0;
+// }

ZZ_test/error_plot.py

@@ -0,0 +1,35 @@
+import numpy as np
+import pandas as pd
+from pylab import *
+import seaborn as sns
+sns.set()
+sns.set_context("paper", font_scale=1.5)
+sns.set_style("ticks")
+
+order = 2
+
+plt.rc('text', usetex=True)
+plt.rc('font', family='serif')
+plt.rc('text.latex', preamble=r'\usepackage{amsmath} \usepackage{amssymb}')
+
+df = pd.read_csv('conv_order%i.csv' % order, index_col=0)
+dofs = df.index.to_numpy()
+H10_errors = df.iloc[:, 0].to_numpy()
+ZZ_errors = df.iloc[:, 1].to_numpy()
+
+H10_error_rate = np.log(H10_errors[-1]/H10_errors[-2])/np.log(dofs[-2]/dofs[-1])
+ZZ_error_rate = np.log(ZZ_errors[-1]/ZZ_errors[-2])/np.log(dofs[-2]/dofs[-1])
+print('H10 error rate: ', H10_error_rate)
+print('Estimator error rate:  ', ZZ_error_rate)
+
+fig1,(ax1) = plt.subplots(1, 1)
+ax1.loglog(dofs,H10_errors,'-og',lw=1.5,markersize = 9.0, alpha=.7, label=r'$H^1_0$ (rate: %.2f)' % H10_error_rate)
+ax1.loglog(dofs,ZZ_errors,'-or',lw=1.5,markersize = 9.0, alpha=.7, label=r'Estimator (rate: %.2f)' % ZZ_error_rate)
+ax1.set_ylabel(r"Error")
+ax1.set_xlabel(r"DOFs")
+ax1.legend(fontsize=15)
+ax1.grid(True)
+
+ax1.set_title(r'Order %i. Expected rate: %.1f' % (order, order/2)  )
+plt.savefig('ConvergencePlotOrder%i.pdf' % order)
+plt.show()

ZZ_test/hp_ZZ_test.cpp

@@ -0,0 +1,496 @@
+//                         Test of ZZ error estimator
+//
+// Compile with: make ZZ_test
+//
+
+#include "mfem.hpp"
+#include <fstream>
+
+// #include "exact.hpp"
+
+using namespace std;
+using namespace mfem;
+
+double lshape_exsol(const Vector &p, double omega);
+void   lshape_exgrad(const Vector &p, double omega, Vector &grad);
+double lshape_laplace(const Vector &p);
+
+int dim;
+const char* keys = "Rjlmc*******";
+
+bool ContainsVertex(Mesh *mesh, int elem, const Vertex& vert) // different name?
+{
+   IsoparametricTransformation Tr;
+   mesh->GetElementTransformation(elem, &Tr);
+   IntegrationPoint reference_pt;
+   Vector physical_pt(3);
+   for (int i = 0; i < 2; i++)
+   {
+      for (int j = 0; j < 2; j++)
+      {
+         reference_pt.Set((float)i, (float)j, 0.0, 0.0);
+         Tr.Transform(reference_pt, physical_pt);
+         double dist = 0.0;
+         for (int l = 0; l < 2; l++)
+         {
+            double d = physical_pt(l) - vert(l);
+            dist += d*d;
+         }
+         if (dist == 0) { return true; }
+      }
+   }
+   return false;
+   // mesh->GetElementVertices(elem, v);
+   // for (int j = 0; j < v.Size(); j++)
+   // {
+   //    double* vertex = mesh->GetVertex(v[j]);
+   //    double dist = 0.0;
+   //    for (int l = 0; l < 2; l++) // Euclidean distance in x-y plane
+   //    {
+   //       double d = vert(l) - vertex[l];
+   //       dist += d*d;
+   //    }
+   //    if (dist == 0) { return true; }
+   // }
+   // return false;
+}
+
+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;
+}
+
+int main(int argc, char *argv[])
+{
+   // Parse command-line options.
+   int problem = 0;
+   int order = 1;
+   double ref_threshold = 0.8;
+   int nc_limit = 1;
+   const char *device_config = "cpu";
+   bool visualization = false;
+   int which_estimator = 0;
+   double angle = 7.0*M_PI/4.0;
+   // double angle = 3.0*M_PI/2.0;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&problem, "-p", "--problem",
+                  "Problem type: 0 = canonical L-shaped solution, 1 = sinusoid.");
+   args.AddOption(&order, "-o", "--order",
+                  "Initial mesh finite element order (polynomial degree).");
+   args.AddOption(&ref_threshold, "-rt", "--ref-threshold",
+                  "Refine elements with error larger than threshold * max_error.");
+   args.AddOption(&angle, "-a", "--angle", "Angle of the reentrant corner.");
+   args.AddOption(&nc_limit, "-nc", "--nc-limit",
+                  "Set maximum difference of refinement levels of adjacent elements.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&which_estimator, "-est", "--estimator",
+                  "Which estimator to use: "
+                  "0 = ZZ, 1 = Kelly. Defaults to ZZ.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   // Enable hardware devices such as GPUs, and programming models such as
+   // CUDA, OCCA, RAJA and OpenMP based on command line options.
+   Device device(device_config);
+   device.Print();
+
+   const char *mesh_file;
+   // if (problem == 0)
+   // {
+      mesh_file = "l-shape-benchmark.mesh";
+   // }
+
+
+   // 2. Read the (serial) mesh from the given mesh file.
+   Mesh mesh(mesh_file, 1, 1);
+   mesh.EnsureNodes();
+   Vector nodes;
+   mesh.GetNodes(nodes);
+   int num_nodes = int(nodes.Size()/2);
+   for (int i = 0; i < num_nodes; i++)
+   {
+      double x = nodes[2*i];
+      double y = nodes[2*i+1];
+      double theta = atan2(y, x);
+      if (theta < 0) { theta += 2.0*M_PI; }
+      double delta_theta = theta * (angle - 3.0*M_PI/2.0) / (3.0*M_PI/2.0);
+      nodes[2*i] = x*cos(delta_theta) - y*sin(delta_theta);
+      nodes[2*i+1] = x*sin(delta_theta) + y*cos(delta_theta);
+   }
+   mesh.SetNodes(nodes);
+
+   dim = mesh.Dimension();
+   mesh.EnsureNCMesh();
+   mesh.UniformRefinement(); // ZZ doesn't work properly on the initial L-shaped mesh
+
+   // Define a finite element space on the mesh.
+   H1_FECollection fec(order, dim);
+   L2_FECollection l2fec(order, dim);
+   FiniteElementSpace fespace(&mesh, &fec);
+
+   // Define the solution vector x as a finite element grid function
+   // corresponding to fespace.
+   GridFunction x(&fespace);
+
+   // Define exact solutions
+   FunctionCoefficient *exsol=nullptr;
+   VectorFunctionCoefficient *exgrad=nullptr;
+   FunctionCoefficient *rhs=nullptr;
+   ConstantCoefficient one(1.0);
+
+//    switch (problem)
+//    {
+//       case 1:
+//       {
+//          exsol = new FunctionCoefficient(sinsin_exsol);
+//          exgrad = new VectorFunctionCoefficient(dim, sinsin_exgrad);
+//          rhs = new FunctionCoefficient(sinsin_laplace);
+//          break;
+//       }
+//       case 2:
+//       {
+//          exsol = new FunctionCoefficient(poly_exsol);
+//          exgrad = new VectorFunctionCoefficient(dim, poly_exgrad);
+//          rhs = new FunctionCoefficient(poly_laplace);
+//          break;
+//       }
+//       default:
+//       case 0:
+//       {
+         exsol = new FunctionCoefficient(lshape_exsol);
+         exgrad = new VectorFunctionCoefficient(dim, lshape_exgrad);
+         rhs = new FunctionCoefficient(lshape_laplace);
+         exsol->SetTime(angle);
+         exgrad->SetTime(angle);
+//          break;
+//       }
+//    }
+
+   // Set up the linear form b(.) and the bilinear form a(.,.).
+   LinearForm b(&fespace);
+   BilinearForm a(&fespace);
+
+   b.AddDomainIntegrator(new DomainLFIntegrator(*rhs));
+   a.AddDomainIntegrator(new DiffusionIntegrator(one));
+
+   // All boundary attributes will be used for essential (Dirichlet) BC.
+   MFEM_VERIFY(mesh.bdr_attributes.Size() > 0,
+               "Boundary attributes required in the mesh.");   
+
+   // Connect to GLVis.
+   socketstream sol_sock, ord_sock, dbg_sock[3], err_sock;
+
+   ostringstream file_name;
+   file_name << "conv_order" << order << ".csv";
+   ofstream conv(file_name.str().c_str());
+//    std::ofstream elemerr(elemerr_file);
+
+   if (visualization)
+   {
+      cout << "\n Press enter to advance... " << endl;
+      cin.get();
+   }
+
+   cout << setw(4) << "\nRef." << setw(12) << "DOFs" << setw(21) << "H^1_0 error" << setw(21) << "error estimate" << setw(18) << "H^1_0 rate" << setw(18) << "estimator rate" <<  endl;
+   conv << "DOFs " << ", " << "H^1_0 error" << ", " << "error estimate" << endl;
+
+   double old_num_dofs = 0.0;
+   double old_H10_error = 0.0;
+   double old_ZZ_error = 0.0;
+   const int max_dofs = 20000;
+   for (int it = 0; ; it++)
+   {
+      int num_dofs = fespace.GetTrueVSize();
+
+      // Set Dirichlet boundary values in the GridFunction x.
+      // Determine the list of Dirichlet true DOFs in the linear system.
+      Array<int> ess_bdr(mesh.bdr_attributes.Max());
+      ess_bdr = 1;
+      x = 0.0;
+      x.ProjectBdrCoefficient(*exsol, ess_bdr);
+      Array<int> ess_tdof_list;
+      fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   
+      // Solve for the current mesh:
+      b.Assemble();
+      a.Assemble();
+      a.Finalize();
+      OperatorPtr A;
+      Vector B, X;
+   
+      const int copy_interior = 1;
+      a.FormLinearSystem(ess_tdof_list, x, b, A, X, B, copy_interior);
+      GSSmoother M((SparseMatrix&)(*A));
+      PCG(*A, M, B, X, 0, 2000, 1e-30, 0.0);
+
+      a.RecoverFEMSolution(X, b, x);
+
+      // Calculate the total error in the H^1_0 norm.
+      double H10_error = x.ComputeGradError(exgrad);
+      DiffusionIntegrator di;
+
+      ErrorEstimator* estimator{nullptr};
+      switch (which_estimator)
+      {
+         case 1:
+         {
+            auto flux_fes = new FiniteElementSpace(&mesh, &fec, dim);
+            estimator = new NewZienkiewiczZhuEstimator(di, x, flux_fes);
+            // int flux_order = 4;
+            // estimator = new NewZienkiewiczZhuEstimator(di, x, flux_order);
+            break;
+         }
+
+         case 2:
+         {
+            auto flux_fes = new FiniteElementSpace(&mesh, &l2fec, dim);
+            estimator = new KellyErrorEstimator(di, x, flux_fes);
+            break;
+         }
+
+         default:
+            std::cout << "Unknown estimator. Falling back to ZZ." << std::endl;
+         case 0:
+         {
+            auto flux_fes = new FiniteElementSpace(&mesh, &fec, dim);
+            estimator = new ZienkiewiczZhuEstimator(di, x, flux_fes);
+            break;
+         }
+      }
+      StopWatch chrono;
+      chrono.Clear();
+      chrono.Start();
+      const Vector &zzerr = estimator->GetLocalErrors();
+      chrono.Stop();
+      mfem::out << "get local errors time: " << chrono.RealTime() << endl;
+      // double ZZ_error = zzerr.Norml2();
+      double ZZ_error = estimator->GetTotalError();
+
+      // estimate convergence rate
+      double H10_rate = 0.0;
+      double ZZ_rate = 0.0;
+      if (old_H10_error > 0.0)
+      {
+          H10_rate = log(H10_error/old_H10_error) / log(old_num_dofs/num_dofs);
+          ZZ_rate  = log(ZZ_error/old_ZZ_error)   / log(old_num_dofs/num_dofs);
+      }
+
+      cout << setw(4) << it << setw(12) << num_dofs << setw(21) << H10_error << setw(21) << ZZ_error << setw(18) << H10_rate << setw(18) << ZZ_rate << endl;
+
+      // Send solution by socket to the GLVis server.
+      if (visualization)
+      {
+         const char vishost[] = "localhost";
+         const int  visport   = 19916;
+
+         // Prolong the solution vector onto L2 space of max order (for GLVis)
+         GridFunction *vis_x = ProlongToMaxOrder(&x);
+         socketstream sol_sock(vishost, visport);
+         sol_sock.precision(8);
+         sol_sock << "solution\n" << mesh << *vis_x;
+         sol_sock << "keys ARjlm\n";
+         delete vis_x;
+
+         // Visualize element orders
+         if (true)
+         {
+            L2_FECollection l20fec(0, dim);
+            FiniteElementSpace l20fes(&mesh, &l20fec);
+            GridFunction orders(&l20fes);
+
+            for (int i = 0; i < orders.Size(); i++)
+            {
+               orders(i) = fespace.GetElementOrder(i);
+            }
+
+            socketstream ord_sock(vishost, visport);
+            ord_sock.precision(8);
+            ord_sock << "solution\n" << mesh << orders;
+            ord_sock << "keys ARjlmpc**]]]]]]]]]]\n";
+         }
+      }
+
+      if (num_dofs > max_dofs)
+      {
+         cout << "\n Reached the maximum number of dofs. Stop." << endl;
+         break;
+      }
+
+      // Save dofs and error for convergence plot
+      conv << num_dofs << ", " << H10_error << ", " << ZZ_error << endl;
+
+    //   for (int i = 0; i < mesh.GetNE(); i++)
+    //   {
+    //      elemerr << sqrt(elemError[i]) << ' ';
+    //   }
+    //   elemerr << endl;
+
+
+    //   Array<Refinement> refinements;
+    //   double err_max = zzerr.Max();
+    //   for (int i = 0; i < mesh.GetNE(); i++)
+    //   {
+    //      if (zzerr[i] > ref_threshold * err_max)
+    //      {
+    //          refinements.Append(Refinement(i, 7));
+    //      }
+    //   }
+    //   mesh.GeneralRefinement(refinements, -1, nc_limit);
+      int h_refined = 0, p_refined = 0;
+
+      Array<Refinement> h_refinements;
+      Array<int> p_refinements;
+      const Table& table = mesh.ElementToElementTable();
+      const Vertex origin(0,0);
+
+      double err_max = zzerr.Max();
+      for (int i = 0; i < mesh.GetNE(); i++)
+      {
+         if (zzerr[i] > ref_threshold * err_max)
+         {
+            if (ContainsVertex(&mesh, i, origin))
+            {
+               h_refinements.Append(Refinement(i));
+               h_refined++;
+               cout << "=> h-refined" << endl;
+            }
+            else
+            {
+               p_refinements.Append(i);
+               // Also refine face-neighbors up to this order
+               int elem_p = fespace.GetElementOrder(i);
+               const int* row = table.GetRow(i);
+               int row_size = table.RowSize(i);
+               for (int j = 0; j < row_size; j++)
+               {
+                  int neig_p = fespace.GetElementOrder(row[j]);
+                  if (neig_p <= elem_p)
+                  {
+                     p_refinements.Append(row[j]);
+                  }
+               }
+            }
+         }
+      }
+      // Filter for unique elements
+      p_refinements.Sort();
+      p_refinements.Unique();
+      for (auto i: p_refinements)
+      {
+         int p = fespace.GetElementOrder(i);
+         fespace.SetElementOrder(i, p+1);
+         p_refined++;
+         cout << "=> p-refined" << endl;
+      }
+
+      old_num_dofs = double(num_dofs);
+      old_H10_error = H10_error;
+      old_ZZ_error = ZZ_error;
+
+      // Update the space, interpolate the solution.
+      fespace.Update(false);
+      mesh.GeneralRefinement(h_refinements, -1, nc_limit);
+      fespace.Update(false);
+
+      a.Update();
+      b.Update();
+      x.Update();
+
+      if (visualization)
+      {
+         cin.get();
+      }
+
+      // Free the used memory.
+      delete estimator;
+   }
+
+   // Free the used memory.
+   delete exsol;
+   delete exgrad;
+   delete rhs;
+   return 0;
+}
+
+
+// L-shape domain problem exact solution (2D)
+
+double lshape_exsol(const Vector &p, double omega)
+{
+   double alpha = M_PI / omega;
+   double x = p(0), y = p(1);
+   double r = sqrt(x*x + y*y);
+   double t = atan2(y, x);
+   if (t < 0) { t += 2.0*M_PI; }
+   return pow(r, alpha) * sin(t*alpha);
+}
+
+void lshape_exgrad(const Vector &p, double omega, Vector &grad)
+{
+   double alpha = M_PI / omega;
+   double x = p(0), y = p(1);
+   double r = sqrt(x*x + y*y);
+   double t = atan2(y, x);
+   if (t < 0) { t += 2*M_PI; }
+   double talpha = t*alpha;
+   double ralpha = pow(r, alpha - 2.0);
+   grad(0) = alpha*ralpha*(x*sin(talpha) - y*cos(talpha));
+   grad(1) = alpha*ralpha*(y*sin(talpha) + x*cos(talpha));
+}
+
+double lshape_laplace(const Vector &p)
+{
+   return 0.0;
+}

ZZ_test/l-shape-benchmark.mesh

@@ -0,0 +1,44 @@
+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
+3
+1 3 0 3 4 1
+1 3 3 6 7 4
+1 3 4 5 2 1
+
+boundary
+8
+1 1 0 1
+1 1 1 2
+1 1 2 5
+2 1 5 4
+2 1 4 7
+1 1 7 6
+1 1 6 3
+1 1 3 0
+
+vertices
+8
+2
+-1 1
+0 1
+1 1
+-1 0
+0 0
+1 0
+-1 -1
+0 -1

ZZ_test/makefile

@@ -0,0 +1,110 @@
+
+
+# Use the MFEM build directory
+MFEM_DIR ?= ..
+MFEM_BUILD_DIR ?= ..
+SRC = $(if $(MFEM_DIR:..=),$(MFEM_DIR)/ZZ_test/,)
+CONFIG_MK = $(MFEM_BUILD_DIR)/config/config.mk
+# Use the MFEM install directory
+# MFEM_INSTALL_DIR = ../mfem
+# CONFIG_MK = $(MFEM_INSTALL_DIR)/share/mfem/config.mk
+
+MFEM_LIB_FILE = mfem_is_not_built
+-include $(CONFIG_MK)
+
+SEQ_EXAMPLES = zz_test
+PAR_EXAMPLES = 
+SEQ_DEVICE_EXAMPLES = 
+PAR_DEVICE_EXAMPLES = 
+
+ifeq ($(MFEM_USE_MPI),NO)
+   EXAMPLES = $(SEQ_EXAMPLES)
+else
+   EXAMPLES = $(PAR_EXAMPLES) $(SEQ_EXAMPLES)
+endif
+SUBDIRS =
+ifeq ($(MFEM_USE_AMGX),YES)
+   SUBDIRS += amgx
+endif
+ifeq ($(MFEM_USE_GINKGO),YES)
+   SUBDIRS += ginkgo
+endif
+ifeq ($(MFEM_USE_HIOP),YES)
+   SUBDIRS += hiop
+endif
+ifeq ($(MFEM_USE_PETSC),YES)
+   SUBDIRS += petsc
+endif
+ifeq ($(MFEM_USE_PUMI),YES)
+   SUBDIRS += pumi
+endif
+ifeq ($(MFEM_USE_SUNDIALS),YES)
+   SUBDIRS += sundials
+endif
+ifeq ($(MFEM_USE_SUPERLU),YES)
+   SUBDIRS += superlu
+endif
+ifeq ($(MFEM_USE_CALIPER),YES)
+   SUBDIRS += caliper
+endif
+
+SUBDIRS_ALL = $(addsuffix /all,$(SUBDIRS))
+SUBDIRS_TEST = $(addsuffix /test,$(SUBDIRS))
+SUBDIRS_CLEAN = $(addsuffix /clean,$(SUBDIRS))
+SUBDIRS_TPRINT = $(addsuffix /test-print,$(SUBDIRS))
+
+.SUFFIXES:
+.SUFFIXES: .o .cpp .mk
+.PHONY: all clean clean-build clean-exec
+
+# Remove built-in rule
+%: %.cpp
+
+# Replace the default implicit rule for *.cpp files
+%: $(SRC)%.cpp $(MFEM_LIB_FILE) $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) $< -o $@ $(MFEM_LIBS)
+
+all: $(EXAMPLES) $(SUBDIRS_ALL)
+
+.PHONY: $(SUBDIRS_ALL) $(SUBDIRS_TEST) $(SUBDIRS_CLEAN) $(SUBDIRS_TPRINT)
+$(SUBDIRS_ALL) $(SUBDIRS_TEST) $(SUBDIRS_CLEAN):
+	$(MAKE) -C $(@D) $(@F)
+$(SUBDIRS_TPRINT):
+	@$(MAKE) -C $(@D) $(@F)
+
+MFEM_TESTS = EXAMPLES
+include $(MFEM_TEST_MK)
+test: $(SUBDIRS_TEST)
+test-print: $(SUBDIRS_TPRINT)
+
+%.o: $(SRC)%.cpp $(MFEM_LIB_FILE) $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) -c $< -o $@
+
+# ZZ_test: ZZ_test.o exact.o $(MFEM_LIB_FILE)
+# 	$(MFEM_CXX) $(MFEM_FLAGS) $(<) -o $(@) ZZ_test.o exact.o $(MFEM_LIB_FILE)
+
+# Generate an error message if the MFEM library is not built and exit
+$(MFEM_LIB_FILE):
+	$(error The MFEM library is not built)
+
+clean: clean-build clean-exec $(SUBDIRS_CLEAN)
+
+clean-build:
+	rm -f *.o *~ $(SEQ_EXAMPLES) $(PAR_EXAMPLES)
+	rm -rf *.dSYM *.TVD.*breakpoints
+
+# clean-exec:
+# 	@rm -f refined.mesh displaced.mesh mesh.* ex5.mesh ex6p-checkpoint.*
+# 	@rm -rf Example5* Example9* Example15* Example16* Example23* ParaView
+# 	@rm -f sphere_refined.* sol.* sol_u.* sol_p.* sol_r.* sol_i.*
+# 	@rm -f ex9.mesh ex9-mesh.* ex9-init.* ex9-final.*
+# 	@rm -f deformed.* velocity.* elastic_energy.* mode_* flux.*
+# 	@rm -f ex5-p-*.bp ex9-p-*.bp ex12-p-*.bp ex16-p-*.bp
+# 	@rm -f ex16.mesh ex16-mesh.* ex16-init.* ex16-final.*
+# 	@rm -f vortex-mesh.* vortex.mesh vortex-?-init.* vortex-?-final.*
+# 	@rm -f deformation.* pressure.*
+# 	@rm -f ex20.dat ex20p_?????.dat gnuplot_ex20.inp gnuplot_ex20p.inp
+# 	@rm -f ex21*.mesh ex21*.sol ex21p_*.*
+# 	@rm -f ex23.mesh ex23-*.gf
+# 	@rm -f ex25.mesh ex25-*.gf ex25p-*.*
+# 	@rm -rf ex28_* ex28p_*

examples/ex15.cpp

@@ -109,7 +109,7 @@ int main(int argc, char *argv[])
                   "Final time; start time is 0.");
    args.AddOption(&which_estimator, "-est", "--estimator",
                   "Which estimator to use: "
-                  "0 = ZZ, 1 = Kelly. Defaults to ZZ.");
+                  "0 = ZZ, 1 = Kelly, 2 = True ZZ. Defaults to ZZ.");
    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                   "--no-visualization",
                   "Enable or disable GLVis visualization.");
@@ -213,6 +213,12 @@ int main(int argc, char *argv[])
          estimator = new KellyErrorEstimator(*integ, x, flux_fes);
          break;
       }
+      case 2:
+      {
+         auto flux_fes = new FiniteElementSpace(&mesh, &fec, sdim);
+         estimator = new NewZienkiewiczZhuEstimator(*integ, x, flux_fes);
+         break;
+      }
 
       default:
          std::cout << "Unknown estimator. Falling back to ZZ." << std::endl;

examples/ex6.cpp

@@ -94,6 +94,7 @@ int main(int argc, char *argv[])
    // 4. Since a NURBS mesh can currently only be refined uniformly, we need to
    //    convert it to a piecewise-polynomial curved mesh. First we refine the
    //    NURBS mesh a bit more and then project the curvature to quadratic Nodes.
+   // mesh.UniformRefinement();
    if (mesh.NURBSext)
    {
       for (int i = 0; i < 2; i++)
@@ -152,8 +153,11 @@ int main(int argc, char *argv[])
    //     flux to get an error indicator. We need to supply the space for the
    //     smoothed flux: an (H1)^sdim (i.e., vector-valued) space is used here.
    FiniteElementSpace flux_fespace(&mesh, &fec, sdim);
-   ZienkiewiczZhuEstimator estimator(*integ, x, flux_fespace);
-   estimator.SetAnisotropic();
+   // ZienkiewiczZhuEstimator estimator(*integ, x, flux_fespace);
+
+   NewZienkiewiczZhuEstimator estimator(*integ, x, flux_fespace);
+
+   estimator.SetAnisotropic(false);
 
    // 11. A refiner selects and refines elements based on a refinement strategy.
    //     The strategy here is to refine elements with errors larger than a
@@ -161,6 +165,7 @@ int main(int argc, char *argv[])
    //     The refiner will call the given error estimator.
    ThresholdRefiner refiner(estimator);
    refiner.SetTotalErrorFraction(0.7);
+   // refiner.SetNCLimit(1);
 
    // 12. The main AMR loop. In each iteration we solve the problem on the
    //     current mesh, visualize the solution, and refine the mesh.
@@ -220,6 +225,7 @@ int main(int argc, char *argv[])
       // 19. Send solution by socket to the GLVis server.
       if (visualization && sol_sock.good())
       {
+         cout << "num elements = " << mesh.GetNE() << endl;
          sol_sock.precision(8);
          sol_sock << "solution\n" << mesh << x << flush;
       }
@@ -230,11 +236,23 @@ int main(int argc, char *argv[])
          break;
       }
 
+      // ParaViewDataCollection paraview_dc("Example5", &mesh);
+      // paraview_dc.SetPrefixPath("ParaView");
+      // paraview_dc.SetLevelsOfDetail(order);
+      // paraview_dc.SetCycle(0);
+      // paraview_dc.SetDataFormat(VTKFormat::BINARY);
+      // paraview_dc.SetHighOrderOutput(true);
+      // paraview_dc.SetTime(0.0); // set the time
+      // paraview_dc.RegisterField("solution",&x);
+      // paraview_dc.Save();
+
       // 20. Call the refiner to modify the mesh. The refiner calls the error
       //     estimator to obtain element errors, then it selects elements to be
       //     refined and finally it modifies the mesh. The Stop() method can be
       //     used to determine if a stopping criterion was met.
       refiner.Apply(mesh);
+      cout << "error = " << estimator.GetTotalError() << endl;
+      // cin.get();
       if (refiner.Stop())
       {
          cout << "Stopping criterion satisfied. Stop." << endl;

examples/face_neigbors.cpp

@@ -0,0 +1,75 @@
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   const char *mesh_file = "../data/inline-quad.mesh";
+   bool vis = true;
+   OptionsParser args(argc, argv);
+   args.AddOption(&vis, "-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_file, 1, 1);
+
+   Array<int> ref_elems(1);
+   ref_elems = 2;
+   mesh.GeneralRefinement(ref_elems,-1,0);
+   ref_elems = 5;
+   mesh.GeneralRefinement(ref_elems,-1,0);
+
+   <<<<<<< HEAD
+   Array<int> faces({0,1,10,26,28,29});
+   =======
+      Array<int> faces({0,10,26,28,29});
+   >>>>>>> zz_test
+
+   for (int i = 0; i<faces.Size(); i++)
+   {
+      Array<int> elems;
+      <<<<<<< HEAD
+      int type = mesh.GetFaceElements2(faces[i],elems);
+      cout << "type = " << type << endl;
+      cout << "elems = " ; elems.Print();
+      =======
+         mesh.GetFaceElements(faces[i],elems);
+
+      >>>>>>> zz_test
+      if (vis)
+      {
+         char vishost[] = "localhost";
+         int  visport   = 19916;
+         socketstream patch_sock(vishost, visport);
+         L2_FECollection fec(1,mesh.Dimension());
+         FiniteElementSpace fes(&mesh,&fec);
+         GridFunction vis_gf(&fes);
+         vis_gf = 0.0;
+         Array<int> dofs;
+         for (int i=0; i<elems.Size(); i++)
+         {
+            int el = elems[i];
+            fes.GetElementDofs(el, dofs);
+            vis_gf.SetSubVector(dofs,1.0);
+         }
+         patch_sock.precision(8);
+         patch_sock << "solution\n" << mesh << vis_gf <<
+                    "keys rRmjnppppp\n" << flush;
+      }
+   }
+
+   return 0;
+}

fem/bilininteg.cpp

@@ -964,7 +964,8 @@ void DiffusionIntegrator::AssembleElementVector(
 
 void DiffusionIntegrator::ComputeElementFlux
 ( const FiniteElement &el, ElementTransformation &Trans,
-  Vector &u, const FiniteElement &fluxelem, Vector &flux, bool with_coef )
+  Vector &u, const FiniteElement &fluxelem, Vector &flux, bool with_coef,
+  const IntegrationRule *ir)
 {
    int i, j, nd, dim, spaceDim, fnd;
 
@@ -1001,13 +1002,17 @@ void DiffusionIntegrator::ComputeElementFlux
    vecdxt.SetSize(spaceDim);
    pointflux.SetSize(MQ || VQ ? spaceDim : 0);
 
-   const IntegrationRule &ir = fluxelem.GetNodes();
-   fnd = ir.GetNPoints();
+   // const IntegrationRule &ir = fluxelem.GetNodes();
+   if (!ir)
+   {
+      ir = &fluxelem.GetNodes();
+   }
+   fnd = ir->GetNPoints();
    flux.SetSize( fnd * spaceDim );
 
    for (i = 0; i < fnd; i++)
    {
-      const IntegrationPoint &ip = ir.IntPoint(i);
+      const IntegrationPoint &ip = ir->IntPoint(i);
       el.CalcDShape(ip, dshape);
       dshape.MultTranspose(u, vec);
 
@@ -1942,7 +1947,7 @@ void CurlCurlIntegrator::AssembleElementMatrix
 void CurlCurlIntegrator
 ::ComputeElementFlux(const FiniteElement &el, ElementTransformation &Trans,
                      Vector &u, const FiniteElement &fluxelem, Vector &flux,
-                     bool with_coef)
+                     bool with_coef, const IntegrationRule *ir)
 {
 #ifdef MFEM_THREAD_SAFE
    DenseMatrix projcurl;
@@ -2776,7 +2781,7 @@ void ElasticityIntegrator::AssembleElementMatrix(
 void ElasticityIntegrator::ComputeElementFlux(
    const mfem::FiniteElement &el, ElementTransformation &Trans,
    Vector &u, const mfem::FiniteElement &fluxelem, Vector &flux,
-   bool with_coef)
+   bool with_coef, const IntegrationRule *ir)
 {
    const int dof = el.GetDof();
    const int dim = el.GetDim();
@@ -2799,14 +2804,18 @@ void ElasticityIntegrator::ComputeElementFlux(
    DenseMatrix gh(gh_data, dim, dim);
    DenseMatrix grad(grad_data, dim, dim);
 
-   const IntegrationRule &ir = fluxelem.GetNodes();
-   const int fnd = ir.GetNPoints();
+   // const IntegrationRule &ir = fluxelem.GetNodes();
+   if (!ir)
+   {
+      ir = &fluxelem.GetNodes();
+   }
+   const int fnd = ir->GetNPoints();
    flux.SetSize(fnd * tdim);
 
    DenseMatrix loc_data_mat(u.GetData(), dof, dim);
    for (int i = 0; i < fnd; i++)
    {
-      const IntegrationPoint &ip = ir.IntPoint(i);
+      const IntegrationPoint &ip = ir->IntPoint(i);
       el.CalcDShape(ip, dshape);
       MultAtB(loc_data_mat, dshape, gh);
 

fem/bilininteg.hpp

@@ -213,12 +213,17 @@ public:
                              of the method may choose not to scale the "flux"
                              function by any coefficients describing the
                              integrator.
+       @param[in] ir  If passed (the defualt value is NULL), the implementation
+                        of the method will ignore the fluxelem parameter and,
+                        instead, compute the discrete flux at the points specified
+                        by the integration rule ir.
     */
    virtual void ComputeElementFlux(const FiniteElement &el,
                                    ElementTransformation &Trans,
                                    Vector &u,
                                    const FiniteElement &fluxelem,
-                                   Vector &flux, bool with_coef = true) { }
+                                   Vector &flux, bool with_coef = true,
+                                   const IntegrationRule *ir = NULL) { }
 
    /** @brief Virtual method required for Zienkiewicz-Zhu type error estimators.
 
@@ -2029,7 +2034,8 @@ public:
    virtual void ComputeElementFlux(const FiniteElement &el,
                                    ElementTransformation &Trans,
                                    Vector &u, const FiniteElement &fluxelem,
-                                   Vector &flux, bool with_coef = true);
+                                   Vector &flux, bool with_coef = true,
+                                   const IntegrationRule *ir = NULL);
 
    virtual double ComputeFluxEnergy(const FiniteElement &fluxelem,
                                     ElementTransformation &Trans,
@@ -2448,7 +2454,8 @@ public:
    virtual void ComputeElementFlux(const FiniteElement &el,
                                    ElementTransformation &Trans,
                                    Vector &u, const FiniteElement &fluxelem,
-                                   Vector &flux, bool with_coef);
+                                   Vector &flux, bool with_coef,
+                                   const IntegrationRule *ir = NULL);
 
    virtual double ComputeFluxEnergy(const FiniteElement &fluxelem,
                                     ElementTransformation &Trans,
@@ -2772,7 +2779,8 @@ public:
                                    ElementTransformation &Trans,
                                    Vector &u,
                                    const FiniteElement &fluxelem,
-                                   Vector &flux, bool with_coef = true);
+                                   Vector &flux, bool with_coef = true,
+                                   const IntegrationRule *ir = NULL);
 
    /** Compute the element energy (integral of the strain energy density)
        corresponding to the stress represented by @a flux which is a vector of

fem/estimators.cpp

@@ -30,6 +30,23 @@ void ZienkiewiczZhuEstimator::ComputeEstimates()
    current_sequence = solution->FESpace()->GetMesh()->GetSequence();
 }
 
+void NewZienkiewiczZhuEstimator::ComputeEstimates()
+{
+   flux_space->Update(false);
+   // In parallel, 'flux' can be a GridFunction, as long as 'flux_space' is a
+   // ParFiniteElementSpace and 'solution' is a ParGridFunction.
+   GridFunction flux(flux_space);
+
+   if (!anisotropic) { aniso_flags.SetSize(0); }
+   total_error = NewZZErrorEstimator(*integ, *solution, flux,
+                                     error_estimates,
+                                     flux_averaging,
+                                     with_coeff,
+                                     tichonov_coeff);
+
+   current_sequence = solution->FESpace()->GetMesh()->GetSequence();
+}
+
 
 #ifdef MFEM_USE_MPI
 
@@ -317,23 +334,6 @@ void KellyErrorEstimator::ComputeEstimates()
 
    current_sequence = solution->FESpace()->GetMesh()->GetSequence();
 
-#ifdef MFEM_USE_MPI
-   if (!isParallel)
-#endif // MFEM_USE_MPI
-   {
-      // Finalize element errors
-      for (int e = 0; e < xfes->GetNE(); e++)
-      {
-         auto factor = compute_element_coefficient(mesh, e);
-         // The sqrt belongs to the norm and hₑ to the indicator.
-         error_estimates(e) = sqrt(factor * error_estimates(e));
-      }
-
-      total_error = error_estimates.Norml2();
-      delete flux;
-      return;
-   }
-
 #ifdef MFEM_USE_MPI
 
    // 3. Add error contribution from shared interior faces
@@ -456,7 +456,20 @@ void KellyErrorEstimator::ComputeEstimates()
    MPI_Allreduce(&process_local_error, &total_error, 1, MPI_DOUBLE,
                  MPI_SUM, pfes->GetComm());
    total_error = sqrt(total_error);
+   return;
 #endif // MFEM_USE_MPI
+
+   // Finalize element errors (if serial)
+   for (int e = 0; e < xfes->GetNE(); e++)
+   {
+      auto factor = compute_element_coefficient(mesh, e);
+      // The sqrt belongs to the norm and hₑ to the indicator.
+      error_estimates(e) = sqrt(factor * error_estimates(e));
+   }
+
+   total_error = error_estimates.Norml2();
+   delete flux;
+
 }
 
 void LpErrorEstimator::ComputeEstimates()

fem/estimators.hpp

@@ -166,10 +166,9 @@ public:
 
    /** @brief Set the way the flux is averaged (smoothed) across elements.
 
-       When @a fa is zero (default), averaging is performed globally. When @a fa
-       is non-zero, the flux averaging is performed locally for each mesh
-       attribute, i.e. the flux is not averaged across interfaces between
-       different mesh attributes. */
+       When @a fa is zero (default), averaging is performed across interfaces
+       between different mesh attributes. When @a fa is non-zero, the flux is
+       not averaged across interfaces between different mesh attributes. */
    void SetFluxAveraging(int fa) { flux_averaging = fa; }
 
    /// Return the total error from the last error estimate.
@@ -202,6 +201,172 @@ public:
    }
 };
 
+/** @brief The NewZienkiewiczZhuEstimator class implements the Zienkiewicz-Zhu
+    error estimation procedure [1,2] using face-based patches [3].
+
+    [1] Zienkiewicz, O.C. and Zhu, J.Z., The superconvergent patch recovery
+    and a posteriori error estimates. Part 1: The recovery technique.
+    Int. J. Num. Meth. Engng. 33, 1331-1364 (1992).
+
+    [2] Zienkiewicz, O.C. and Zhu, J.Z., The superconvergent patch recovery
+    and a posteriori error estimates. Part 2: Error estimates and adaptivity.
+    Int. J. Num. Meth. Engng. 33, 1365-1382 (1992).
+
+    [3] Bartels, S. and Carstensen, C., Each averaging technique yields reliable
+    a posteriori error control in FEM on unstructured grids. Part II: Higher
+    order FEM. Math. Comp. 71(239), 971-994 (2002)
+
+    The required BilinearFormIntegrator must implement the methods
+    ComputeElementFlux() and ComputeFluxEnergy().
+
+   COMMENTS:
+   *  The present implementation ignores all single-element patches corresponding
+      to boundary faces. This is appropriate for Dirichlet boundaries, but
+      suboptimal for Neumann boundaries. Reference 3 shows that a constrained
+      least-squares problem, where the reconstructed flux is constrained by the
+      Neumann boundary data, is appropriate to handle this case.
+      THIS CONSTRAINED LS PROBLEM IS NOT YET IMPLEMENTED, so it is possible that
+      the local error estimates for Neumann boundary elements may be affected.
+   *  THIS IMPLEMENTATION IS ONLY SERIAL.
+   *  ANISOTROPIC REFINEMENT NOT YET SUPPORTED.
+
+ */
+class NewZienkiewiczZhuEstimator : public AnisotropicErrorEstimator
+{
+protected:
+   long current_sequence;
+   Vector error_estimates;
+   double total_error;
+   bool anisotropic;
+   Array<int> aniso_flags;
+   int flux_averaging; // see SetFluxAveraging()
+   double tichonov_coeff;
+
+   BilinearFormIntegrator *integ; ///< Not owned.
+   GridFunction *solution; ///< Not owned.
+
+   FiniteElementSpace *flux_space; /**< @brief Ownership based on own_flux_fes.
+      Its Update() method is called automatically by this class when needed. */
+   bool with_coeff;
+   bool own_flux_fes; ///< Ownership flag for flux_space.
+
+   /// Check if the mesh of the solution was modified.
+   bool MeshIsModified()
+   {
+      long mesh_sequence = solution->FESpace()->GetMesh()->GetSequence();
+      MFEM_ASSERT(mesh_sequence >= current_sequence, "");
+      return (mesh_sequence > current_sequence);
+   }
+
+   /// Compute the element error estimates.
+   void ComputeEstimates();
+
+public:
+   /** @brief Construct a new NewZienkiewiczZhuEstimator object.
+    * The arguments are intentionally similar to the ZienkiewiczZhuEstimator
+    * constructor
+       @param integ    This BilinearFormIntegrator must implement the methods
+                       ComputeElementFlux() and ComputeFluxEnergy().
+       @param sol      The solution field whose error is to be estimated.
+       @param flux_fes The ZienkiewiczZhuEstimator assumes ownership of this
+                       FiniteElementSpace and will call its Update() method when
+                       needed.*/
+   NewZienkiewiczZhuEstimator(BilinearFormIntegrator &integ, GridFunction &sol,
+                              FiniteElementSpace *flux_fes)
+      : current_sequence(-1),
+        total_error(),
+        anisotropic(false),
+        flux_averaging(0),
+        tichonov_coeff(0.0),
+        integ(&integ),
+        solution(&sol),
+        flux_space(flux_fes),
+        with_coeff(false),
+        own_flux_fes(true)
+   { }
+
+   /** @brief Construct a new NewZienkiewiczZhuEstimator object.
+    * The arguments are intentionally similar to the ZienkiewiczZhuEstimator
+    * constructor
+       @param integ    This BilinearFormIntegrator must implement the methods
+                       ComputeElementFlux() and ComputeFluxEnergy().
+       @param sol      The solution field whose error is to be estimated.
+       @param flux_fes The ZienkiewiczZhuEstimator does NOT assume ownership of
+                       this FiniteElementSpace; will call its Update() method
+                       when needed. */
+   NewZienkiewiczZhuEstimator(BilinearFormIntegrator &integ, GridFunction &sol,
+                              FiniteElementSpace &flux_fes)
+      : current_sequence(-1),
+        total_error(),
+        anisotropic(false),
+        flux_averaging(0),
+        tichonov_coeff(0.0),
+        integ(&integ),
+        solution(&sol),
+        flux_space(&flux_fes),
+        with_coeff(false),
+        own_flux_fes(false)
+   { }
+
+   /** @brief Consider the coefficient in BilinearFormIntegrator to calculate
+       the fluxes for the error estimator.*/
+   void SetWithCoeff(bool w_coeff = true) { with_coeff = w_coeff; }
+
+   /** @brief Enable/disable anisotropic estimates. To enable this option, the
+       BilinearFormIntegrator must support the 'd_energy' parameter in its
+       ComputeFluxEnergy() method. */
+   void SetAnisotropic(bool aniso = true)
+   {
+      MFEM_WARNING("Anisotropic refinement is not implemented yet.")
+      anisotropic = aniso;
+   }
+
+   /** @brief Solve a Tichonov-regularized least-squares problem for the
+    *         reconstructed fluxes. This is epsecially helpful for when not
+    *         using tensor product elements, which typically require fewer
+    *         integration points. */
+   void SetTichonovRegularization(double lambda = 1.0e-8)
+   {
+      tichonov_coeff = lambda;
+   }
+
+   /** @brief Set the way the flux is averaged (smoothed) across elements.
+
+       When @a fa is zero (default), averaging is performed across interfaces
+       between different mesh attributes. When @a fa is non-zero, the flux is
+       not averaged across interfaces between different mesh attributes. */
+   void SetFluxAveraging(int fa) { flux_averaging = fa; }
+
+   /// Return the total error from the last error estimate.
+   virtual double GetTotalError() const override { return total_error; }
+
+   /// Get a Vector with all element errors.
+   virtual const Vector &GetLocalErrors() override
+   {
+      if (MeshIsModified()) { ComputeEstimates(); }
+      return error_estimates;
+   }
+
+   /** @brief Get an Array<int> with anisotropic flags for all mesh elements.
+       Return an empty array when anisotropic estimates are not available or
+       enabled. */
+   virtual const Array<int> &GetAnisotropicFlags() override
+   {
+      if (MeshIsModified()) { ComputeEstimates(); }
+      return aniso_flags;
+   }
+
+   /// Reset the error estimator.
+   virtual void Reset() override { current_sequence = -1; }
+
+   /** @brief Destroy a ZienkiewiczZhuEstimator object. Destroys, if owned, the
+       FiniteElementSpace, flux_space. */
+   virtual ~NewZienkiewiczZhuEstimator()
+   {
+      if (own_flux_fes) { delete flux_space; }
+   }
+};
+
 
 #ifdef MFEM_USE_MPI
 

fem/gridfunc.hpp

@@ -372,6 +372,10 @@ public:
        in each element (not L2 projection). */
    virtual void ProjectCoefficient(Coefficient &coeff);
 
+   /** @brief Project @a coeff Coefficient to @a this GridFunction, only onto
+    * a subset of elements */
+   void ProjectElemCoefficient(Coefficient &coeff, const Array<int> & elems);
+
    /** @brief Project @a coeff Coefficient to @a this GridFunction, using one
        element for each degree of freedom in @a dofs and nodal interpolation on
        that element. */
@@ -486,10 +490,18 @@ public:
                                  const IntegrationRule *irs[] = NULL,
                                  Array<int> *elems = NULL) const;
 
+   /// Returns ||grad u_ex - grad u_h||_L2 in element i for H1 or L2 elements
+   virtual double ComputeElementGradError(int i, VectorCoefficient *exgrad,
+                                          const IntegrationRule *irs[] = NULL) const;
+
    /// Returns ||grad u_ex - grad u_h||_L2 for H1 or L2 elements
    virtual double ComputeGradError(VectorCoefficient *exgrad,
                                    const IntegrationRule *irs[] = NULL) const;
 
+   virtual double ComputeElementGradErrors(VectorCoefficient *exgrad,
+                                           Vector & errors, const Array<int> * elems = nullptr,
+                                           const IntegrationRule *irs[] = NULL) const;
+
    /// Returns ||curl u_ex - curl u_h||_L2 for ND elements
    virtual double ComputeCurlError(VectorCoefficient *excurl,
                                    const IntegrationRule *irs[] = NULL) const;
@@ -529,6 +541,10 @@ public:
    virtual double ComputeH1Error(Coefficient *exsol, VectorCoefficient *exgrad,
                                  const IntegrationRule *irs[] = NULL) const;
 
+   virtual double ComputeElementH1Errors(Coefficient *exsol,
+                                         VectorCoefficient *exgrad, Vector & errors, const Array<int> * elems = nullptr,
+                                         const IntegrationRule *irs[] = NULL) const;
+
    /// Returns the error measured in H(div)-norm for RT elements
    virtual double ComputeHDivError(VectorCoefficient *exsol,
                                    Coefficient *exdiv,
@@ -574,11 +590,18 @@ public:
        the Vector @a error. The result should be of length number of elements,
        for example an L2 GridFunction of order zero using map type VALUE. */
    virtual void ComputeElementLpErrors(const double p, Coefficient &exsol,
-                                       Vector &error,
+                                       Vector &errors,
                                        Coefficient *weight = NULL,
                                        const IntegrationRule *irs[] = NULL
                                       ) const;
 
+   /** Compute the Lp error in one specific element of the mesh. */
+   virtual double ComputeElementLpError(int ielem,
+                                        const double p, Coefficient &exsol,
+                                        Coefficient *weight = NULL,
+                                        const IntegrationRule *irs[] = NULL
+                                       ) const;
+
    virtual void ComputeElementL1Errors(Coefficient &exsol,
                                        Vector &error,
                                        const IntegrationRule *irs[] = NULL
@@ -591,6 +614,11 @@ public:
                                       ) const
    { ComputeElementLpErrors(2.0, exsol, error, NULL, irs); }
 
+
+   virtual void ComputeElementL2Errors(Coefficient &exsol,
+                                       Vector &errors, const Array<int> & elems,
+                                       const IntegrationRule *irs[] = NULL) const;
+
    virtual void ComputeElementMaxErrors(Coefficient &exsol,
                                         Vector &error,
                                         const IntegrationRule *irs[] = NULL
@@ -916,6 +944,83 @@ double ZZErrorEstimator(BilinearFormIntegrator &blfi,
                         int with_subdomains = 1,
                         bool with_coeff = false);
 
+/// Defines the global polynomial space used by NewZZErorrEstimator
+Vector LegendreND(const Vector & x, const Vector &xmax, const Vector &xmin,
+                  int order, int dim, double angle=0.0, const Vector *center=NULL);
+
+/// Defines the a bounding box for the face patches used by NewZZErorrEstimator
+void BoundingBox(Array<int> patch,                // input
+                 FiniteElementSpace *ufes, // input
+                 int order,                // input
+                 Vector &xmin,             // output
+                 Vector &xmax,             // output
+                 double &angle,            // output
+                 Vector &center,           // output
+                 int iface=-1);            // input (optional)
+
+
+
+class PatchLeastSquaresCoefficient : public Coefficient
+{
+private:
+   int dim;
+   Mesh * mesh = nullptr;
+   FiniteElementSpace * fes = nullptr;
+   GridFunction * u = nullptr;
+   int order;
+   Array<int> elems;
+   Vector xmax, xmin;
+   Vector coefficients;
+   void Setup();
+
+public:
+   PatchLeastSquaresCoefficient(GridFunction * u_, int order_,
+                                const Array<int> & elems_)
+      : Coefficient(), u(u_), order(order_), elems(elems_)
+   {
+      fes = u->FESpace();
+      mesh = fes->GetMesh();
+      dim = mesh->Dimension();
+      Setup();
+   }
+
+   double Eval(ElementTransformation &T, const IntegrationPoint &ip)
+   {
+      double x[3];
+      Vector transip(x, 3);
+      T.Transform(ip, transip);
+      Vector p(coefficients.Size());
+      p = LegendreND(transip, xmax, xmin, order, dim);
+      return coefficients*p;
+   }
+};
+
+
+
+/// A ``true'' ZZ error estimator which uses face-based patches
+double NewZZErrorEstimator(BilinearFormIntegrator &blfi,
+                           GridFunction &u,
+                           GridFunction &flux,
+                           Vector &error_estimates,
+                           int with_subdomains = 1,
+                           bool with_coeff = false,
+                           double tichonov_coeff = 0.0);
+
+
+class PatchBasedPolynomialFit
+{
+private:
+   int patch_order, integ_order, dim;
+   Vector coefficients, xmin, xmax;
+public:
+   PatchBasedPolynomialFit(GridFunction &u, Array<int> elems,
+                           int patch_order_, int integ_order_,
+                           double tichonov_coeff);
+
+   double EvaluatePolynomial(const Vector &xloc);
+
+};
+
 /// Compute the Lp distance between two grid functions on the given element.
 double ComputeElementLpDistance(double p, int i,
                                 GridFunction& gf1, GridFunction& gf2);
@@ -939,6 +1044,7 @@ public:
 GridFunction *Extrude1DGridFunction(Mesh *mesh, Mesh *mesh2d,
                                     GridFunction *sol, const int ny);
 
+GridFunction* ProlongToMaxOrder(const GridFunction *x);
 
 // Inline methods
 

mesh/mesh.cpp

@@ -1138,6 +1138,87 @@ void Mesh::GetFaceInfos(int Face, int *Inf1, int *Inf2, int *NCFace) const
    *NCFace = faces_info[Face].NCFace;
 }
 
+void Mesh::GetFaceElements(int face, Array<int> & elems) const
+{
+   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];
+      if (!nc_info.Slave)
+      {
+         const mfem::NCMesh::NCList &nc_list = ncmesh->GetNCList(Dim-1);
+         elems.Append(ncmesh->elements[nc_list.masters[nc_index].element].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);
+         }
+         return;
+      }
+   }
+   // Conforming face or slave face. In case of nonconforming master face, early
+   // return above.
+   int el1, el2;
+   GetFaceElements(face, &el1, &el2);
+   if (el1 != -1) { elems.Append(el1); }
+   if (el2 != -1) { elems.Append(el2); }
+}
+
+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;
+   }
+}
+
 Geometry::Type Mesh::GetFaceGeometryType(int Face) const
 {
    switch (Dim)

mesh/mesh.hpp

@@ -1193,6 +1193,15 @@ public:
    void GetFaceInfos (int Face, int *Inf1, int *Inf2) const;
    void GetFaceInfos (int Face, int *Inf1, int *Inf2, int *NCFace) const;
 
+   /** Return all elements adjacent to the given Face
+       For an NCMesh and a master face the function will return more than 2 elements.
+       For the other two cases (conforming or slave face) it falls back to the
+       original GetFaceElements function above */
+   void GetFaceElements (int Face, Array<int> & elems) const;
+
+   int GetFaceElementsAndFaces(int face, Array<int> & elems,
+                               Array<int> & faces) const;
+
    Geometry::Type GetFaceGeometryType(int Face) const;
    Element::Type  GetFaceElementType(int Face) const;