testing amr results for l-shape domain problem

examples/fd/element-smoother.cpp

@@ -0,0 +1,353 @@
+#include "element-smoother.hpp"
+
+using namespace std;
+using namespace mfem;
+
+
+ElementSmoother::ElementSmoother(ParFiniteElementSpace * fes_,  
+Array<int> ess_bdr, Coefficient * cf_ )
+: Solver(fes_->GetTrueVSize()), fes(fes_), cf(cf_)
+{  
+   comm = fes->GetComm();
+   MPI_Comm_size(comm, &num_procs);
+   MPI_Comm_rank(comm, &myid);
+   Pr = fes->GetProlongationMatrix();
+   ParMesh * pmesh = fes->GetParMesh();
+   dim = pmesh->Dimension();
+   eidx.SetSize(dim);
+   eidx[0] = 0;
+   eidx[1] = 1;
+   if (dim == 3) eidx[2] = 8;
+   nrelems = pmesh->GetNE();
+   tpcf.SetSize(nrelems);
+   int vsize = fes->GetVSize();
+   ovlp_count.SetSize(vsize);
+   ovlp_count = 0.0;
+   // Construct overlap count for each dof &
+   // count x y z edges sharing a vertex
+   for (int i=0; i<nrelems; i++)
+   {
+      Array<int> elem_dofs;
+      fes->GetElementDofs(i,elem_dofs);
+      for (int j = 0; j<elem_dofs.Size(); j++)
+      {
+         ovlp_count(elem_dofs[j]) +=1.0;
+      }
+   }
+
+
+   Vector tovlp_count(fes->TrueVSize());
+   if (Pr)
+   {
+      Pr->MultTranspose(ovlp_count,tovlp_count);
+      Pr->Mult(tovlp_count, ovlp_count);
+   }
+
+   DenseMatrix edge_counts;
+   // GetVertexToEdgeCount(pmesh,edge_counts);
+   double * data = edge_counts.GetData();
+
+   // helper H1 fespace for communication of vertex info
+   H1_FECollection fec(1, dim); 
+   ParFiniteElementSpace aux_fes(const_cast<ParMesh *>(pmesh), &fec);
+   Vector tedge_counts(aux_fes.TrueVSize());
+   Vector temp(aux_fes.GetVSize());
+
+
+   // const Operator * Ph = aux_fes.GetProlongationMatrix();
+
+   // if (Ph)
+   // {
+   //    for (int d=0; d<dim; d++)
+   //    {
+   //       temp.SetData(&data[d*aux_fes.GetVSize()]);
+   //       Ph->MultTranspose(temp,tedge_counts);
+   //       Ph->Mult(tedge_counts, temp);
+   //    }
+   // }
+
+   for (int i = 0; i < vsize; i++)
+   {
+      ovlp_count(i) = 1.0/sqrt(ovlp_count(i)); 
+   }
+
+   fes->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   TPElementTransformation TPTrans(*fes);
+   
+   int nredges = pmesh->GetNEdges();
+   edge_orient.SetSize(nredges); // orientation of edges wrt the element
+   Array<DenseMatrix *> EdgeGrad(nredges);
+   Array<DenseMatrix *> EdgeMass(nredges);
+   Array<DenseMatrix *> AssembledEdgeGrad(nredges);
+   Array<DenseMatrix *> AssembledEdgeMass(nredges);
+   // Initialize
+   for (int i = 0; i<nredges; i++)
+   {
+      AssembledEdgeGrad[i] = nullptr;
+      AssembledEdgeMass[i] = nullptr;
+      EdgeGrad[i] = nullptr;
+      EdgeMass[i] = nullptr;
+   }
+
+   // loop through element to calculate Edge matrices
+   Array<int> emarker(nredges); emarker = 0;
+   for (int iel = 0; iel<nrelems; iel++)
+   {
+      tpcf[iel] = new ElementTPFunctionCoefficient(*fes,iel,*cf);
+      Array<int> edges, cor;
+      pmesh->GetElementEdges(iel,edges,cor);
+      for (int ii = 0; ii<dim; ii++) 
+      {   
+         int i = eidx[ii];
+         int edge = edges[i];
+         if (emarker[edge]) continue;
+         edge_orient[edge] = cor[i];
+         const FiniteElement * fe = fes->GetEdgeElement(edge);
+
+         tpcf[iel]->SetCoord(ii);
+         tpcf[iel]->SetOrient(edge_orient[edge]);
+         tpcf[iel]->ResetCounter(ii);
+         IntegrationRule *irs = TensorIntegrationRule(1,fe->GetOrder());
+         EdgeGrad[edge] = new DenseMatrix(fe->GetDof());
+         EdgeMass[edge] = new DenseMatrix(fe->GetDof());
+
+         int j;
+         Vector * Q;
+         switch (ii)
+         {
+         case 0: 
+            j = 1; 
+            Q = tpcf[iel]->GetVecX();
+            break;
+         case 1: 
+            j = 0; 
+            Q = tpcf[iel]->GetVecY();
+            break;  
+         default:
+            j=2; 
+            Q = tpcf[iel]->GetVecZ();
+            break;
+         }
+
+         Vector * vecG = TPTrans.GetTPTransformation(iel,ii,ii);
+         Vector * vecM = nullptr;
+         Vector *vecM1 = nullptr;
+         Vector *vecM2 = nullptr;
+         if (dim == 2)
+         {
+            vecM = TPTrans.GetTPTransformation(iel,ii,j);
+         }
+         else
+         {
+            switch (ii)
+            {
+            case 0:
+            {
+               vecM1 = TPTrans.GetTPTransformation(iel,ii,1);
+               vecM2 = TPTrans.GetTPTransformation(iel,ii,2);
+            }
+            break;
+            case 1:
+            {
+               vecM1 = TPTrans.GetTPTransformation(iel,ii,0);
+               vecM2 = TPTrans.GetTPTransformation(iel,ii,2);
+            }            
+            default:
+            {
+               vecM1 = TPTrans.GetTPTransformation(iel,ii,0);
+               vecM2 = TPTrans.GetTPTransformation(iel,ii,1);
+            }
+               break;
+            }
+            vecM = new Vector(vecM1->Size());
+            for (int i=0; i<vecM->Size(); i++)
+            {
+               // (*vecM)(i) = ((*vecM1)(i)+(*vecM2)(i))/2.0;
+               // (*vecM)(i) = (*vecM1)(i);
+               (*vecM)(i) = 1.0/((1./(*vecM1)(i)+1./(*vecM2)(i))/2.0);
+            }   
+         }
+         // GetDiffusionEdgeMatrix(edge,fes,*vecG,*Q,irs,*EdgeGrad[edge],edge_orient[edge]);
+         // GetMassEdgeMatrix(edge,fes,*vecM,*Q,irs,*EdgeMass[edge],edge_orient[edge]);
+         Get1DMatrices(fes,edge, edge_orient[edge],
+                       *vecG, *vecM,*Q,irs,*EdgeGrad[edge],*EdgeMass[edge]);
+         emarker[edge] = 1;
+      }
+   }
+
+   
+   for (int iel=0; iel<nrelems; iel++)
+   {
+      Array<int> edges, cor;
+      pmesh->GetElementEdges(iel,edges,cor);
+
+      for (int ii = 0; ii<dim; ii++) 
+      {   
+         int i = eidx[ii];
+         int k = edges[i];
+         const FiniteElement *fe = fes->GetEdgeElement(k);
+         int ndof = fe->GetDof();
+         DenseMatrix Grad(ndof); Grad = *EdgeGrad[k];
+         DenseMatrix Mass(ndof); Mass = *EdgeMass[k];
+
+         Array<int> vert;
+         pmesh->GetEdgeVertices(k,vert);
+         for (int i = 0; i<2; i++)
+         {
+            int vertex = vert[i];
+            Array<int> vertdofs;
+            fes->GetVertexDofs(vertex,vertdofs);
+            // double count = edge_counts(vertdofs[0],ii); 
+
+            // Grad(i,i) *= count;
+            // Mass(i,i) *= count;
+            Grad(i,i) *= 2.;
+            Mass(i,i) *= 2.;
+         }
+         const Array<int> &dmap =
+         dynamic_cast<const TensorBasisElement&>(*fe).GetDofMap();
+         SparseMatrix * P = new SparseMatrix(dmap.Size());
+         for (int j = 0; j<dmap.Size(); j++)
+         {
+            P->Set(dmap[j],j, 1.0);
+         }
+         P->Finalize();
+         // Map from MFEM ordering to TensorProduct Ordering
+         AssembledEdgeGrad[k] = RAP(Grad,*P);
+         AssembledEdgeMass[k] = RAP(Mass,*P);
+         delete P;
+      }
+   }
+
+   Array<Array<int> * > tmap;
+   TensorProductEssentialDofsMaps(ess_tdof_list, fes, tmap, dofmap);
+
+
+   Array<DenseMatrix * > G(nredges);
+   Array<DenseMatrix * > M(nredges);
+
+   for (int ie = 0; ie<nredges; ie++)
+   {
+      if (!emarker[ie]) continue;
+      if (!AssembledEdgeMass[ie])
+      {
+         cout << "ie = " << ie << endl;
+         MFEM_ABORT("Memory allocation incosistency 2");
+      }
+      const FiniteElement * fe = fes->GetEdgeElement(ie);
+      int n = fe->GetDof() - tmap[ie]->Size();
+      G[ie] = new DenseMatrix(n);
+      M[ie] = new DenseMatrix(n);
+      const Array<int> &dmap =
+      dynamic_cast<const TensorBasisElement&>(*fe).GetDofMap();
+      // modify tmap to use tensor product index;
+      Array<int> dmapt(dmap.Size());
+      for (int i = 0; i<dmap.Size(); i++)
+      {
+         dmapt[dmap[i]] = i;
+      }
+
+      // Eliminate indices corresponding to tmap from matrices Grad1D_A and Mass1D_A
+      // construct Map;
+      Array<int> tmap_marker(fe->GetDof());
+      tmap_marker = 0;
+      for (int i = 0; i<tmap[ie]->Size(); i++)
+      {
+         int j = (*tmap[ie])[i];
+         tmap_marker[j] = 1;
+      }
+      Array<int> dof_list;
+      for (int i =0; i<fe->GetDof(); i++)
+      {
+         if (tmap_marker[i]) continue;
+         dof_list.Append(dmapt[i]);
+      }
+      dof_list.Sort();
+      for (int i=0; i<dof_list.Size(); i++)
+      {
+         int iii = (edge_orient[ie] == 1) ? i : n - i - 1; 
+         int ii = dof_list[i];
+         for (int j=0; j<dof_list.Size(); j++)
+         {
+            int jjj = (edge_orient[ie] == 1) ? j : n - j - 1; 
+            int jj = dof_list[j];
+            (*G[ie])(iii,jjj) = (*AssembledEdgeGrad[ie])(ii,jj);
+            (*M[ie])(iii,jjj) = (*AssembledEdgeMass[ie])(ii,jj);
+         }
+      }
+   }
+
+   elem_inv.SetSize(nrelems);
+   for (int iel = 0; iel<nrelems; iel++)
+   {
+      elem_inv[iel] = nullptr;
+      Array<int> edges, cor;
+      pmesh->GetElementEdges(iel,edges,cor);
+      Array<DenseMatrix *> Gv(dim);
+      Array<DenseMatrix *> Mv(dim);
+      bool msize = true;
+      for (int d = 0; d<dim; d++)
+      {
+         Gv[dim - d - 1] = G[edges[eidx[d]]];
+         Mv[dim - d - 1] = M[edges[eidx[d]]];
+         if (Mv[dim-d-1]->Size() == 0) msize = false;
+      }
+      if (msize) elem_inv[iel] = new FDSolver(Gv,Mv);
+   }
+}
+
+
+void ElementSmoother::Mult(const Vector &r, Vector &z) const
+{
+   z = r;
+   z.SetSubVectorComplement(ess_tdof_list,0.0);
+   Vector rnew(fes->GetVSize());
+   Vector znew(fes->GetVSize()); 
+   Vector ztemp(fes->GetTrueVSize()); 
+   ztemp = 0.0;
+   znew = 0.0;
+
+   // const SparseMatrix * R = fes->GetRestrictionMatrix();
+   if (Pr)
+   {
+      Pr->Mult(r,rnew);
+   }
+   else
+   {
+      rnew = r;
+   }
+   for (int iel=0; iel<nrelems; iel++)
+   {
+      if (!elem_inv[iel]) continue;
+      int n = dofmap[iel]->Size();
+
+      Vector rloc(n);
+      rnew.GetSubVector(*dofmap[iel],rloc);
+      // pre-scale
+      for (int i = 0; i<n; i++)
+      {
+         int j = (*dofmap[iel])[i];
+         rloc[i] *= ovlp_count[j]; 
+      }
+      Vector zloc(n);
+      elem_inv[iel]->Mult(rloc,zloc);
+      // post-scale
+      for (int i = 0; i<n; i++)
+      {
+         int j = (*dofmap[iel])[i];
+         zloc[i] *= ovlp_count[j]; 
+      }
+
+      znew.AddElementVector(*dofmap[iel],zloc);
+   }
+   if (Pr)
+   {
+      Pr->MultTranspose(znew,ztemp);
+   }
+   else
+   {
+      ztemp = znew;
+   }
+   // R->Mult(znew,ztemp);
+   z += ztemp;
+}

examples/fd/element-smoother.hpp

@@ -0,0 +1,30 @@
+#pragma once
+#include "smoother-util.hpp"
+
+using namespace std;
+using namespace mfem;
+
+class ElementSmoother: public Solver
+{
+private:
+   int num_procs, myid;
+   MPI_Comm comm;
+   int nrelems;
+   int dim;
+   ParFiniteElementSpace * fes = nullptr;
+   const Operator * Pr = nullptr;
+   Coefficient * cf = nullptr;
+   Array<int> eidx; // edge local index
+   Array<int> edge_orient; // orientation of edges wrt the element
+   Array<FDSolver *> elem_inv;
+   Array<int> ess_tdof_list;
+   Vector ovlp_count;
+   Array<Array<int> * > dofmap;
+   Array<ElementTPFunctionCoefficient *> tpcf;
+public:
+   ElementSmoother(ParFiniteElementSpace * fes_, Array<int> ess_bdr, Coefficient * cf_=nullptr);
+   virtual void SetOperator(const Operator &op) { }
+   virtual void Mult(const Vector &r, Vector &z) const;
+   virtual void MultTranspose(const Vector &r, Vector &z) const { Mult(r,z); }
+   virtual ~ElementSmoother(){};
+};

examples/fd/ex_amr_diffusionp.cpp

@@ -0,0 +1,289 @@
+//                                MFEM Example 26
+//
+
+#include "mfem.hpp"
+#include "exact_sol.hpp"
+#include <fstream>
+#include <iostream>
+#include <random>
+#include "element-smoother.hpp"
+using namespace std;
+using namespace mfem;
+
+int main(int argc, char *argv[])
+{
+   // 0. Initialize MPI.
+   int num_procs, myid;
+   MPI_Init(&argc, &argv);
+   MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
+   MPI_Comm_rank(MPI_COMM_WORLD, &myid);
+
+   // 1. Parse command-line options.
+   const char *mesh_file = "l-shape-benchmark.mesh";
+   int init_geometric_refinements = 0;
+   int pinit_geometric_refinements = 0;
+   int geometric_refinements = 0;
+   int order_refinements = 2;
+   const char *device_config = "cpu";
+   bool visualization = true;
+   int order = 1;
+   int solver = 0;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order", "Finite element order.");
+   args.AddOption(&solver, "-solver", "--solver", "Solver: 0:MG-Cheb-Jac, 1: MG-Cheb-ElemSmoother");
+   args.AddOption(&init_geometric_refinements, "-ref", "--initial-geometric-refinements",
+                  "Number of serial geometric refinements defining the coarse mesh.");
+   args.AddOption(&pinit_geometric_refinements, "-pref", "--initial-geometric-refinements",
+                  "Number of parallel geometric refinements defining the coarse mesh.");                  
+   args.AddOption(&geometric_refinements, "-gr", "--geometric-refinements",
+                  "Number of geometric refinements done prior to order refinements.");
+   args.AddOption(&order_refinements, "-or", "--order-refinements",
+                  "Number of order refinements. Finest level in the hierarchy has order 2^{or}.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      MPI_Finalize();
+      return 1;
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+
+   Mesh *mesh = new Mesh(mesh_file, 1, 1);
+   int dim = mesh->Dimension();
+
+   for (int l = 0; l < init_geometric_refinements; l++)
+   {
+      mesh->UniformRefinement();
+   }
+   mesh->EnsureNCMesh();
+   ParMesh * pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
+   mesh->Clear();
+   {
+      for (int l = 0; l < pinit_geometric_refinements; l++)
+      {
+         pmesh->UniformRefinement();
+      }
+   }
+
+   FiniteElementCollection *fec = new H1_FECollection(order, dim);
+   ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, fec);
+
+   Array<int> ess_bdr;
+   if(pmesh->bdr_attributes.Size())
+   {   
+      ess_bdr.SetSize(pmesh->bdr_attributes.Max());
+      ess_bdr = 1;
+   }
+
+
+   ParGridFunction x(fespace);
+   FunctionCoefficient ex_coeff(lshape_exact);
+   x.ProjectCoefficient(ex_coeff);
+
+   // ------------------------------------------------- 
+   // Bilinear and linear forms
+   // ------------------------------------------------- 
+   ConstantCoefficient cf(1.0);
+   ParBilinearForm a(fespace);
+   a.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+   DiffusionIntegrator * aa = new DiffusionIntegrator(cf);
+   // int order1 = fespace->GetElementOrder(0);
+   IntegrationRule *irs = TensorIntegrationRule(*fespace,order); 
+   aa->SetIntegrationRule(*irs);
+   a.AddDomainIntegrator(aa);
+
+   ParLinearForm b(fespace);
+   FunctionCoefficient rhscf(lshape_rhs);
+   b.AddDomainIntegrator(new DomainLFIntegrator(rhscf));
+   // ------------------------------------------------- 
+
+   char vishost[] = "localhost";
+   int  visport   = 19916;
+
+   socketstream sout;
+   if (visualization)
+   {
+      sout.open(vishost, visport);
+      if (!sout)
+      {
+         if (myid == 0)
+         {
+            cout << "Unable to connect to GLVis server at "
+                 << vishost << ':' << visport << endl;
+            cout << "GLVis visualization disabled.\n";
+         }
+         visualization = false;
+      }
+
+      sout.precision(8);
+
+      sout << "parallel " << num_procs << " " << myid << "\n";
+      sout << "solution\n" << *pmesh << x << flush;
+   }
+
+   L2_FECollection flux_fec(order, dim);
+   ParFiniteElementSpace flux_fes(pmesh, &flux_fec, dim);
+   FiniteElementCollection *smooth_flux_fec = NULL;
+   ParFiniteElementSpace *smooth_flux_fes = NULL;
+   smooth_flux_fec = new RT_FECollection(order-1, dim);
+   smooth_flux_fes = new ParFiniteElementSpace(pmesh, smooth_flux_fec, 1);
+   L2ZienkiewiczZhuEstimator estimator(*aa, x, flux_fes, *smooth_flux_fes);
+   ThresholdRefiner refiner(estimator);
+   refiner.SetTotalErrorFraction(0.7);
+   refiner.SetNCLimit(1);
+   StopWatch chrono;
+   Array<double> ts0, ts1, tsol;
+   int ref_amr = 20;   
+   Array<int> iter;
+   Array<int> dofs;
+
+   ostringstream file_name;
+   file_name << "lshape-amr_" << order << ".csv";
+   ofstream conv(file_name.str().c_str());
+   conv << "DOFs " << ", " << "it-Cheb-Jac" << ", " << "it-ChebElemSmoother" << endl;
+
+   for (int it = 0; it < ref_amr ; it++)
+   {
+      HYPRE_BigInt global_dofs = fespace->GlobalTrueVSize();
+      if (myid == 0)
+      {
+         cout << "\nAMR iteration " << it << endl;
+         cout << "Number of unknowns: " << global_dofs << endl;
+      }
+
+      Array<int> ess_tdof_list;
+      fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+
+      a.Assemble();
+      b.Assemble();
+
+      OperatorPtr A;
+      Vector B, X;
+      a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
+
+
+      chrono.Clear();
+      chrono.Start();
+      chrono.Stop();
+      ts0.Append(chrono.RealTime());
+      chrono.Clear();
+      chrono.Start();
+      Solver * prec1 = nullptr;
+      Solver * prec2 = nullptr;
+      Solver * S = nullptr;
+      // if (solver)
+      // {
+         S = new ElementSmoother(fespace,ess_bdr, &cf); 
+         prec1 = new OperatorChebyshevSmoother(*A, *S, 4, MPI_COMM_WORLD,7);
+      // }
+      // else
+      // {
+         Vector diag(fespace->GetTrueVSize());
+         a.AssembleDiagonal(diag);
+         prec2 = new OperatorChebyshevSmoother(*A, diag,ess_tdof_list, 4, MPI_COMM_WORLD,10);
+      // }
+
+      chrono.Stop();
+      ts1.Append(chrono.RealTime());
+
+
+
+      int print_level = 3;
+      int max_iter = 2000;
+      double rtol = 1e-8;
+
+      CGSolver pcg(MPI_COMM_WORLD);
+      pcg.SetPrintLevel(print_level);
+      pcg.SetMaxIter(max_iter);
+      pcg.SetRelTol(rtol);
+      pcg.SetOperator(*A);
+      pcg.SetPreconditioner(*prec1);
+
+
+      // chrono.Clear();
+      // chrono.Start();
+      Vector Y = X;
+      pcg.Mult(B,Y);
+      int iter1 = pcg.GetNumIterations();
+
+
+      pcg.SetPreconditioner(*prec2);
+      pcg.Mult(B,X);
+      int iter2 = pcg.GetNumIterations();
+
+      // chrono.Stop();
+      // tsol.Append(chrono.RealTime());
+      // iter.Append(pcg.GetNumIterations());
+      // dofs.Append(global_dofs);
+
+      delete S; 
+      delete prec1;
+      delete prec2;
+
+      conv << global_dofs << ", " << iter1 << ", " << iter2 << endl;
+
+
+      a.RecoverFEMSolution(X,b,x);
+
+      if (visualization)
+      {
+         sout << "parallel " << num_procs << " " << myid << "\n";
+         sout << "solution\n" << *pmesh << x << flush;
+      }
+      refiner.Apply(*pmesh);
+      if (refiner.Stop())
+      {
+         if (myid == 0)
+         {
+            cout << "Stopping criterion satisfied. Stop." << endl;
+         }
+         break;
+      }
+      fespace->Update();
+      x.Update();
+      x.ProjectCoefficient(ex_coeff);
+
+      a.Update();
+      b.Update();
+   }
+
+
+
+   if (myid==0)
+   {
+      cout << "ts0 total  = " << ts0.Sum() << endl;
+      cout << "ts1 total  = " << ts1.Sum() << endl;
+      cout << "tsol total = " << tsol.Sum() << endl;
+   }
+
+   if (myid == 0)
+   {
+      cout << "num iterations = "; iter.Print(cout, iter.Size());
+      cout << "dofs = "; dofs.Print(cout, dofs.Size());
+   }
+
+   delete smooth_flux_fes;
+   delete smooth_flux_fec;
+   delete pmesh;
+
+   // 13. Free the used memory.
+   MPI_Finalize();
+   return 0;
+}
+
+
+

examples/fd/ex_diffusionp.cpp

@@ -0,0 +1,670 @@
+//                                MFEM Example 26
+//
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+#include <random>
+#include "element-smoother.hpp"
+using namespace std;
+using namespace mfem;
+
+class DiffusionMultigrid : public GeometricMultigrid
+{
+private:
+   Coefficient * cf = nullptr;
+   int smoother_kind = 0;
+   // 0: Jacobi, 1:Chebychev, 2: element-smoother(matrix-free)
+   HypreBoomerAMG* amg;
+
+public:
+   // Constructs a diffusion multigrid for the ParFiniteElementSpaceHierarchy
+   // and the array of essential boundaries
+   DiffusionMultigrid(ParFiniteElementSpaceHierarchy& fespaces,
+                      Array<int>& ess_bdr, Coefficient * cf_,int smoother_ = 0)
+      : GeometricMultigrid(fespaces), cf(cf_), smoother_kind(smoother_)
+   {
+      ConstructCoarseOperatorAndSolver(fespaces.GetFESpaceAtLevel(0), ess_bdr,cf);
+
+      for (int level = 1; level < fespaces.GetNumLevels(); ++level)
+      {
+         ConstructOperatorAndSmoother(fespaces.GetFESpaceAtLevel(level), ess_bdr,cf);
+      }
+   }
+
+   virtual ~DiffusionMultigrid()
+   {
+      delete amg;
+   }
+
+private:
+   void ConstructBilinearForm(ParFiniteElementSpace& fespace, Array<int>& ess_bdr,
+                              bool partial_assembly, Coefficient * cf)
+   {
+      ParBilinearForm* form = new ParBilinearForm(&fespace);
+      if (partial_assembly)
+      {
+         form->SetAssemblyLevel(AssemblyLevel::PARTIAL);
+      }
+      form->AddDomainIntegrator(new DiffusionIntegrator(*cf));
+      form->Assemble();
+      bfs.Append(form);
+
+      essentialTrueDofs.Append(new Array<int>());
+      fespace.GetEssentialTrueDofs(ess_bdr, *essentialTrueDofs.Last());
+   }
+
+   void ConstructCoarseOperatorAndSolver(ParFiniteElementSpace& coarse_fespace,
+                                         Array<int>& ess_bdr, 
+                                         Coefficient * cf)
+   {
+      ConstructBilinearForm(coarse_fespace, ess_bdr, false, cf);
+
+      HypreParMatrix* hypreCoarseMat = new HypreParMatrix();
+      bfs.Last()->FormSystemMatrix(*essentialTrueDofs.Last(), *hypreCoarseMat);
+
+      amg = new HypreBoomerAMG(*hypreCoarseMat);
+      amg->SetPrintLevel(-1);
+
+      CGSolver* pcg = new CGSolver(MPI_COMM_WORLD);
+      pcg->SetPrintLevel(-1);
+      pcg->SetMaxIter(10);
+      pcg->SetRelTol(sqrt(1e-8));
+      pcg->SetAbsTol(0.0);
+      pcg->SetOperator(*hypreCoarseMat);
+      pcg->SetPreconditioner(*amg);
+
+      AddLevel(hypreCoarseMat, pcg, true, true);
+   }
+
+   void ConstructOperatorAndSmoother(ParFiniteElementSpace& fespace,
+                                     Array<int>& ess_bdr, Coefficient *cf)
+   {
+      ConstructBilinearForm(fespace, ess_bdr, true, cf);
+
+      OperatorPtr opr;
+      opr.SetType(Operator::ANY_TYPE);
+      bfs.Last()->FormSystemMatrix(*essentialTrueDofs.Last(), opr);
+      opr.SetOperatorOwner(false);
+      Solver * smoother = nullptr;
+      Vector diag;
+      if (smoother_kind < 2 )
+      {
+         diag.SetSize(fespace.GetTrueVSize());
+         bfs.Last()->AssembleDiagonal(diag);
+      }
+
+      switch (smoother_kind)
+      {
+      case 0:
+         smoother = new OperatorJacobiSmoother(diag,*essentialTrueDofs.Last(),0.6667);
+         break;
+      case 1:
+         smoother = new OperatorChebyshevSmoother(opr.Ptr(), diag,
+                                       *essentialTrueDofs.Last(), 5);
+         break;   
+      case 2:
+         smoother = new ElementSmoother(&fespace,ess_bdr,cf);
+         break;      
+      case 3:
+         {
+            ElementSmoother * sm = new ElementSmoother(&fespace,ess_bdr,cf);
+            smoother = new OperatorChebyshevSmoother(*opr,*sm,5,fespace.GetComm());
+            break;      
+         }
+      default:
+         MFEM_ABORT("Wrong smoother choice");
+         break;
+      }
+      AddLevel(opr.Ptr(), smoother, true, true);
+   }
+};
+
+int dim;
+int  exact   = 0;
+bool tpcoeff = true;
+
+double f_exact(const Vector & x);
+double u_exact(const Vector & x);
+void usol(const Vector & x, double &u, Vector & Grad, double & d2u);
+
+double DiffusionCoeff(const Vector & x);
+double TPDiffusionCoeff(const Vector & x, int coord);
+void DiffusionCoeffGrad(const Vector & x, Vector & Grad);
+
+
+int main(int argc, char *argv[])
+{
+   // 0. Initialize MPI.
+   int num_procs, myid;
+   MPI_Init(&argc, &argv);
+   MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
+   MPI_Comm_rank(MPI_COMM_WORLD, &myid);
+
+   // 1. Parse command-line options.
+   const char *mesh_file = "../../data/inline-quad.mesh";
+   int init_geometric_refinements = 0;
+   int pinit_geometric_refinements = 0;
+   int geometric_refinements = 0;
+   int order_refinements = 2;
+   const char *device_config = "cpu";
+   bool visualization = true;
+   int order = 1;
+   double skew_factor = 0.0;
+   double scale_factor = 1.0;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order", "Finite element order.");
+   args.AddOption(&init_geometric_refinements, "-ref", "--initial-geometric-refinements",
+                  "Number of serial geometric refinements defining the coarse mesh.");
+   args.AddOption(&pinit_geometric_refinements, "-pref", "--initial-geometric-refinements",
+                  "Number of parallel geometric refinements defining the coarse mesh.");                  
+   args.AddOption(&geometric_refinements, "-gr", "--geometric-refinements",
+                  "Number of geometric refinements done prior to order refinements.");
+   args.AddOption(&order_refinements, "-or", "--order-refinements",
+                  "Number of order refinements. Finest level in the hierarchy has order 2^{or}.");
+   args.AddOption(&tpcoeff, "-tpcoeff", "--tp-coefficient", "-no-tpcoeff",
+               "--no-tp-coefficient", "Tensor product diffusion coefficient or not");
+   args.AddOption(&exact, "-exact", "--exact", "Exact Solution flag: 0: unknown");   
+   args.AddOption(&skew_factor, "-c", "--skew_factor", "Skew_factor");   
+   args.AddOption(&scale_factor, "-s", "--scale_factor", "Scale_factor");   
+   
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      MPI_Finalize();
+      return 1;
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+
+   // 2. 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);
+   if (myid == 0) { device.Print(); }
+
+
+   // 3. Read the mesh from the given mesh file. We can handle triangular,
+   //    quadrilateral, tetrahedral, hexahedral, surface and volume meshes with
+   //    the same code.
+   Mesh *mesh = new Mesh(mesh_file, 1, 1);
+   // int nx = pow(2,init_geometric_refinements);
+   // int ny = pow(2,init_geometric_refinements);
+   // Mesh *mesh = new Mesh(1,1,1,mfem::Element::HEXAHEDRON,true,1.0,2.0,3.0,false);
+   // Mesh *mesh = new Mesh(1,1,mfem::Element::QUADRILATERAL,true,1.0,1.0,false);
+   // Mesh *mesh = new Mesh(1,4.0);
+   // move nodes
+   dim = mesh->Dimension();
+   mesh->EnsureNodes();
+   
+   mesh->SetCurvature(3);
+   GridFunction * nodes = mesh->GetNodes();
+   // *nodes +=1.0;
+   // *nodes *=0.5;
+   double c = skew_factor;
+   double s = scale_factor;
+   if (dim == 2)
+   {
+      for (int i=0; i<nodes->Size()/2; i++)
+      {
+      // double temp = (*nodes)(2*i);
+      // (*nodes)(2*i) += (*nodes)(2*i)*(*nodes)(2*i) + c*pow((*nodes)(2*i+1),2);
+      // (*nodes)(2*i) +=  c*pow((*nodes)(2*i+1),2);
+         (*nodes)(2*i) +=  c*pow((*nodes)(2*i+1),2);
+      // (*nodes)(2*i+1) = (*nodes)(2*i+1)*(*nodes)(2*i+1) + c*pow(temp,2);
+      }
+   }
+   else
+   {
+      for (int i=0; i<nodes->Size()/3; i++)
+      {
+         // (*nodes)(3*i) +=  c*pow((*nodes)(3*i+1),2);
+         // (*nodes)(3*i+1) +=  c*pow((*nodes)(3*i+2),3);
+         (*nodes)(3*i+2) +=  c*pow((*nodes)(3*i),2);
+      }
+   }
+
+   for (int i=0; i<nodes->Size(); i++)
+   {
+      (*nodes)(i) *=  s;
+   }
+
+   dim = mesh->Dimension();
+   // mesh->EnsureNCMesh();
+
+   // 4. Refine the mesh to increase the resolution and order
+   {
+      for (int l = 0; l < init_geometric_refinements; l++)
+      {
+         mesh->UniformRefinement();
+      }
+   }
+
+   ParMesh * pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
+   mesh->Clear();
+   {
+      for (int l = 0; l < pinit_geometric_refinements; l++)
+      {
+         pmesh->UniformRefinement();
+      }
+   }
+
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream mesh_sock(vishost, visport);
+      mesh_sock << "parallel " << num_procs << " " << myid << "\n";
+      mesh_sock.precision(8);
+      mesh_sock << "mesh\n" << *pmesh << flush;
+   }
+
+   FiniteElementCollection *fec = new H1_FECollection(order, dim);
+   ParFiniteElementSpace *coarse_fespace = new ParFiniteElementSpace(pmesh, fec);
+   ParFiniteElementSpaceHierarchy fespaces(pmesh, coarse_fespace, true, true);
+
+
+   Coefficient * cf = nullptr;
+   if (exact)
+   {
+      cf = new FunctionCoefficient(DiffusionCoeff);
+   }
+   else
+   {
+      cf = new ConstantCoefficient(1.0);
+   }
+
+   Array<FiniteElementCollection*> collections;
+   collections.Append(fec);
+   for (int level = 0; level < geometric_refinements; ++level)
+   {
+      fespaces.AddUniformlyRefinedLevel();
+   }
+   for (int level = 0; level < order_refinements; ++level)
+   {
+      // order++;
+      order *=2;
+      // collections.Append(new H1_FECollection(std::pow(2, level+1), dim));
+      collections.Append(new H1_FECollection(order, dim));
+      fespaces.AddOrderRefinedLevel(collections.Last());
+   }
+
+   HYPRE_Int size = fespaces.GetFinestFESpace().GlobalTrueVSize();
+
+   if (myid == 0)
+   {
+      cout << "Number of finite element unknowns: " << size << endl;
+      cout << "Order = " << order << endl;
+   }
+
+   // 6. Set up the linear form b(.) which corresponds to the right-hand side of
+   //    the FEM linear system, which in this case is (1,phi_i) where phi_i are
+   //    the basis functions in the finite element fespace.
+   FunctionCoefficient f(f_exact);
+   ConstantCoefficient one(1.0);
+     
+   
+   Array<int> ess_bdr;
+   if(pmesh->bdr_attributes.Size())
+   {   
+      ess_bdr.SetSize(pmesh->bdr_attributes.Max());
+      ess_bdr = 1;
+   }
+   ParGridFunction x(&fespaces.GetFinestFESpace());
+   // GridFunction gf_coeff(&fespaces.GetFinestFESpace());
+   ParMesh * ref_mesh = fespaces.GetFinestFESpace().GetParMesh();
+   L2_FECollection * l2fec = new L2_FECollection(order,dim);
+   ParFiniteElementSpace * l2fes = new ParFiniteElementSpace(ref_mesh,l2fec);
+   ParGridFunction gf_coeff(l2fes);
+   // gf_coeff.ProjectCoefficient(*cf);
+   gf_coeff.ProjectDiscCoefficient(*cf,mfem::GridFunction::AvgType::ARITHMETIC);
+   
+   int print_level = 3;
+   int max_iter = 2000;
+   double rtol = 1e-8;
+   StopWatch chrono;
+   // for (int i = 0; i<=6; i++)
+   for (int i = 0; i<=6; i++)
+   {
+      OperatorPtr A;
+      Vector B, X;
+      Solver * prec = nullptr;
+      CGSolver pcg(MPI_COMM_WORLD);
+      pcg.SetPrintLevel(print_level);
+      pcg.SetMaxIter(max_iter);
+      pcg.SetRelTol(rtol);
+      // i=1; Chebychev-Jacobi-MG
+      // i=2; Chebychev-Element-MG
+      // i=3; Chebychev-Jacobi-Smoother
+      // i=4; Element-Smoother 
+      // i=5; Chebychev-Element-Smoother 
+      ParLinearForm *b = new ParLinearForm(&fespaces.GetFinestFESpace());
+      if (exact)
+      {
+         b->AddDomainIntegrator(new DomainLFIntegrator(f));
+      }
+      else
+      {
+         b->AddDomainIntegrator(new DomainLFIntegrator(one));
+      }
+      b->Assemble();
+      FunctionCoefficient u_ex(u_exact);
+      x = 0.0;
+      if (exact) x.ProjectCoefficient(u_ex);
+
+      if (i<4)
+      {
+         prec = new DiffusionMultigrid(fespaces, ess_bdr, cf,i);
+         dynamic_cast<DiffusionMultigrid *>(prec)->
+                  SetCycleType(Multigrid::CycleType::VCYCLE, 1, 1);
+         dynamic_cast<DiffusionMultigrid *>(prec)->
+                  FormFineLinearSystem(x, *b, A, X, B);
+         if (i == 0)
+         {
+            if (myid == 0)
+               cout << "\nJacobi-MG " << endl;
+         }         
+         else if (i == 1)
+         {
+            if (myid == 0)
+               cout << "\nJacobi-Chebychev-MG " << endl;
+         }
+         else if (i == 2)
+         {
+            if (myid == 0)
+               cout << "\nElement-Smoother-MG " << endl;
+         }
+         else 
+         {
+            if (myid == 0)
+               cout << "\nElement-Chebychev-MG " << endl;
+         }
+         pcg.SetOperator(*A);
+         if (prec) { pcg.SetPreconditioner(*prec); }
+         chrono.Clear();
+         chrono.Start();         
+         pcg.Mult(B,X);
+         chrono.Stop();      
+         if (myid == 0)
+            cout<< "PCG::mult time = " << chrono.RealTime() << endl;       
+        // Recover the solution as a finite element grid function.
+        dynamic_cast<DiffusionMultigrid *>(prec)->RecoverFineFEMSolution(X, *b, x);
+        delete prec;
+      }
+      else
+      {
+         ParBilinearForm a(&fespaces.GetFinestFESpace());
+         a.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+         DiffusionIntegrator * aa = new DiffusionIntegrator(*cf);
+         int order1 = fespaces.GetFinestFESpace().GetOrder(0);
+         IntegrationRule *irs = TensorIntegrationRule(fespaces.GetFinestFESpace(),order1); 
+         aa->SetIntegrationRule(*irs);
+         a.AddDomainIntegrator(aa);
+         a.Assemble();
+         
+         Array<int> ess_tdof_list;
+         fespaces.GetFinestFESpace().GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+         a.FormLinearSystem(ess_tdof_list, x, *b, A, X, B);
+         if (i==4)
+         {
+            Vector diag(fespaces.GetFinestFESpace().GetTrueVSize());
+            a.AssembleDiagonal(diag);
+            prec = new OperatorChebyshevSmoother(A.Ptr(), diag,ess_tdof_list, 1, MPI_COMM_WORLD);
+            if (myid == 0)
+               cout << "\nJacobi-Chebychev " << endl;
+         }
+         else if (i==5)
+         {
+            prec = new ElementSmoother(&fespaces.GetFinestFESpace(),ess_bdr, cf); 
+            if (myid == 0)
+               cout << "\nElementSmoother " << endl;
+         }
+         else
+         {
+            ElementSmoother *S = new ElementSmoother(&fespaces.GetFinestFESpace(),ess_bdr, cf); 
+            prec = new OperatorChebyshevSmoother(*A, *S, 1, MPI_COMM_WORLD);
+            if (myid == 0)
+               cout << "\nElement-Chebychev " << endl;
+         }
+         pcg.SetOperator(*A);
+         if (prec) { pcg.SetPreconditioner(*prec); }
+         chrono.Clear();
+         chrono.Start();         
+         pcg.Mult(B,X);
+         chrono.Stop();
+         if (myid == 0)
+            cout<< "PCG::mult time = " << chrono.RealTime() << endl;        
+         delete prec;
+
+         a.RecoverFEMSolution(X,*b,x);
+      }
+      delete b;
+   }
+
+   // 12. Send the solution by socket to a GLVis server.
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream sol_sock(vishost, visport);
+      sol_sock << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock.precision(8);
+      sol_sock << "solution\n" << *fespaces.GetFinestFESpace().GetMesh() << x <<
+               flush;
+      socketstream coeff_sock(vishost, visport);
+      coeff_sock << "parallel " << num_procs << " " << myid << "\n";
+      coeff_sock.precision(8);
+      coeff_sock << "solution\n" << *fespaces.GetFinestFESpace().GetMesh() << gf_coeff <<
+               flush;         
+   }
+
+   // 13. Free the used memory.
+   for (int level = 0; level < collections.Size(); ++level)
+   {
+      delete collections[level];
+   }
+
+   MPI_Finalize();
+   return 0;
+}
+
+
+double f_exact(const Vector & x)
+{
+   // -div (f * grad(u)) = (f * gradu(0))_x + (f*gradu(1))_y + + (f*gradu(2))_z
+   // = f_x*gradu(0) + f * gradu(0)_x + f_y * gradu(1) + f* gradu(1)_y + f_z * gradu(2) + f* gradu(2)_z
+   // = f_x*gradu(0) + f_y * gradu(1) + f_z * gradu(2) + f*d2u
+   double u;
+   double d2u;
+   Vector gradu;
+   usol(x,u,gradu,d2u);
+   Vector gradf;
+   DiffusionCoeffGrad(x,gradf);
+   double f = DiffusionCoeff(x);
+   double val = gradf * gradu + f*d2u;
+   return -val;
+}
+
+double u_exact(const Vector & x)
+{
+   double u;
+   Vector gradu;
+   double d2u;
+   usol(x,u,gradu,d2u);
+   return u;
+}
+
+void usol(const Vector & x, double &u, Vector & Grad, double & d2u)
+{
+   Grad.SetSize(dim);
+   if (exact == 1)
+   {
+      Vector alpha(dim); alpha = 5.0;
+      // Vector alpha(dim); alpha = 0.5;
+      double s = alpha * x; // dot product
+      u = sin(M_PI*s);
+      d2u = 0.0;
+      for (int i = 0; i<dim; i++)
+      {
+         Grad[i] = alpha(i) * M_PI * cos(M_PI*s);
+         d2u += alpha(i)*alpha(i); 
+      }
+      d2u = - M_PI*M_PI * d2u * u;
+   }
+   else if (exact == 2)
+   {
+      double c_0 = 1.2;
+      double k_0 = 3.0;
+      double c_1 = 2.3;
+      double k_1 = 5.0;
+      double c_2 = 1.3;
+      double k_2 = 1.0;
+
+      double alpha = c_0 + k_0 * x(0);
+      double beta  = c_1 + k_1 * x(1);
+      double gamma = 1.0;
+      
+      if (dim == 2)
+      {
+         u = sin(M_PI * alpha) *  sin(M_PI * beta);
+         Grad[0] = M_PI*k_0 * cos(alpha) * sin(M_PI*beta);
+         Grad[1] = M_PI*k_1 * cos(beta) * sin(M_PI*alpha);
+      }
+      else if (dim == 3)
+      {
+         gamma  = c_2 + k_2 * x(2);
+         u = sin(M_PI * alpha) *  sin(M_PI * beta) * sin(M_PI*gamma);
+         Grad[0] = M_PI*k_0 * cos(alpha) * sin(M_PI*beta) * sin(M_PI*gamma);
+         Grad[1] = M_PI*k_1 * cos(beta) * sin(M_PI*alpha) * sin(M_PI*gamma);
+         Grad[2] = M_PI*k_2 * cos(gamma) * sin(M_PI*alpha) * sin(M_PI*beta);
+      }
+
+      double u_xx = - M_PI * M_PI * k_0 * k_0 * u;
+      double u_yy = - M_PI * M_PI * k_1 * k_1 * u;
+      double u_zz = - M_PI * M_PI * k_2 * k_2 * u;
+      d2u = u_xx + u_yy; 
+      if (dim == 3 ) d2u += u_zz;
+   }
+}
+
+double TPDiffusionCoeff(const Vector & x, int coord)
+{
+   double val;
+   switch (coord)
+   {
+      case  0: val = 4.+3.*x(0); break;
+      case  1: val = 0.5+7.*x(1)*x(1); break;
+      case  2: val = (0.1+2.*x(2)); break;
+      default: 
+         val = (4.+3.*x(0))*(0.5+7.*x(1)*x(1)); 
+         if (dim == 3 ) val *= (0.1+2.*x(2));
+      break;
+      // case  0: val = x(0); break;
+      // case  1: val = 1.0; break;
+      // case  2: val = 1.0; break;
+      // default: val = x(0); break;
+      // case  0: val = 3.0; break;
+      // case  1: val = 2.0; break;
+      // case  2: val = 1.0; break;
+      // default: val = 6.0; break;
+   }
+   return val;
+   // return 2.0;
+
+}
+
+double DiffusionCoeff(const Vector & x)
+{
+   double val;
+   if (tpcoeff)
+   {
+      val = (4.+3.*x(0))*(0.5+7.*x(1)*x(1));
+      if (dim == 3) val *= (0.1+2.*x(2));
+   }
+   else
+   {
+      // val = 2.0+cos(x.Sum());
+      Vector cf(dim);
+      // cf(0) = 0.1; cf(1) = 3.; 
+      cf(0) = 1.0; cf(1) = 2.0; 
+      // if (dim == 3) cf(2) = -7.8;
+      if (dim == 3) cf(2) = +1.8;
+      // double dd = x * cf + 1.5* x(1)*x(1);
+      double dd = x * cf;
+      // // double dd = x * cf;
+      // // val = exp(cos(dd));
+      val = exp(dd);
+
+      // Vector alpha(dim); alpha = 5.0;
+      // double s = alpha * x; // dot product
+      // val = 2.0+sin(M_PI*s);
+
+   }   
+   return val;
+}
+
+
+void DiffusionCoeffGrad(const Vector & x, Vector & Grad)
+{
+   Grad.SetSize(dim);
+   if (tpcoeff)
+   {
+      if (dim == 2)
+      {
+         Grad[0] = 3.* (0.5+7.*x(1)*x(1));
+         Grad[1] = 14.* x(1) * (4.+3.*x(0));
+      }
+      else
+      {
+         Grad[0] = 3.* (0.5+7.*x(1)*x(1))*(0.1+2.*x(2));
+         Grad[1] = 14.* x(1) * (4.+3.*x(0))*(0.1+2.*x(2));
+         Grad[2] = 2.*(4.+3.*x(0))*(0.5+7.*x(1)*x(1));
+      }
+   }
+   else
+   {
+      Vector cf(dim);
+      // cf(0) = 0.1; cf(1) = 3.; 
+      cf(0) = 1.0; cf(1) = 2.0; 
+      // if (dim == 3) cf(2) = -7.8;
+      if (dim == 3) cf(2) = 1.8;
+      // double dd = x * cf  + 1.5* x(1)*x(1);
+      double dd = x * cf;
+      Vector alpha(dim); alpha = 5.0;
+
+      // for (int d = 0; d<dim; d++)
+      // {
+      // //    // Grad[d] = -sin(x.Sum());
+      // //    // Grad[d] = -cf(d) * exp(cos(dd))*sin(dd);
+      // //    Grad[d] = cf(d) * exp(dd);
+      //    Grad[d] = alpha(d) * M_PI * cos(M_PI*s);
+      // }
+      if (dim == 2)
+      {
+         Grad[0] = (cf(0) )*exp(dd);
+         Grad[1] = (cf(1) + 3.0*x(1))*exp(dd);
+      }
+      else
+      {
+         // Grad[0] = (cf(0) + 1.5 * x(1))*exp(dd);
+         Grad[0] = cf(0)*exp(dd);
+         Grad[1] = cf(1)*exp(dd);
+         Grad[2] = cf(2)*exp(dd);
+      }
+   }
+}
+
+

examples/fd/exact_sol.cpp

@@ -0,0 +1,39 @@
+
+#include "exact_sol.hpp"
+
+double lshape_exact(const Vector & pt)
+{
+   double x = pt[0];
+   double y = pt[1];
+   double r = sqrt(x*x + y*y);
+   double alpha = 2. / 3.;
+
+   double theta = atan2(y, x);
+   if (y < 0) { theta += 2 * M_PI; }
+
+   return pow(r,alpha) * sin(alpha * theta);
+}
+
+void lshape_grad(const Vector & x, Vector & grad)
+{
+
+}
+
+double lshape_rhs(const Vector & x)
+{
+   return 0.0;
+}
+
+double wavefront_exact(const Vector & x)
+{
+   return 0.;
+}
+
+void wavefront_grad(const Vector & x, Vector & grad)
+{
+
+}
+double wavefront_rhs(const Vector & x)
+{
+   return 0.;
+}

examples/fd/exact_sol.hpp

@@ -0,0 +1,15 @@
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+double lshape_exact(const Vector & x);
+void lshape_grad(const Vector & x, Vector & grad);
+double lshape_rhs(const Vector & x);
+
+double wavefront_exact(const Vector & x);
+void wavefront_grad(const Vector & x, Vector & grad);
+double wavefront_rhs(const Vector & x);

examples/fd/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

examples/fd/makefile

@@ -0,0 +1,59 @@
+# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+# LICENSE and NOTICE for details. LLNL-CODE-806117.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability visit https://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the BSD-3 license. We welcome feedback and contributions, see file
+# CONTRIBUTING.md for details.
+
+
+# Use the MFEM build directory
+MFEM_DIR ?= ../..
+MFEM_BUILD_DIR ?= ../..
+SRC = $(.,)
+CONFIG_MK = $(MFEM_BUILD_DIR)/config/config.mk
+
+MFEM_LIB_FILE = mfem_is_not_built
+-include $(CONFIG_MK)
+
+SEQ_EXAMPLES = ex_diffusion 
+PAR_EXAMPLES = ex_diffusionp  ex_amr_diffusionp
+ifeq ($(MFEM_USE_MPI),NO)
+   EXAMPLES = $(SEQ_EXAMPLES)
+else
+   EXAMPLES = $(PAR_EXAMPLES) $(SEQ_EXAMPLES)
+endif
+
+.SUFFIXES:
+.SUFFIXES: .o .cpp .mk
+.PHONY: all clean
+.PRECIOUS: %.o
+
+COMMON_O= element-smoother.o smoother-util.o exact_sol.o
+
+# Remove built-in rules
+%: %.cpp
+%.o: %.cpp
+
+all: $(EXAMPLES)
+
+# Rules for building the EXAMPLES
+
+%: $(SRC)%.cpp $(COMMON_O) $(MFEM_LIB_FILE) $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) $< -o $@ $(COMMON_O) $(MFEM_LIBS)
+
+# Rules for compiling miniapp dependencies
+$(COMMON_O) $($(EXAMPLES)): \
+%.o: $(SRC)%.cpp $(SRC)%.hpp $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) -c $(<) -o $(@)
+
+# Generate an error message if the MFEM library is not built and exit
+$(MFEM_LIB_FILE):
+	$(error The MFEM library is not built)
+
+clean:
+	rm -f *.o *~ $(SEQ_EXAMPLES) $(PAR_EXAMPLES)
+	rm -rf *.dSYM *.TVD.*breakpoints

examples/fd/smoother-util.cpp

@@ -0,0 +1,746 @@
+
+
+#include "smoother-util.hpp"
+
+IntegrationRule * TensorIntegrationRule(const FiniteElementSpace & fes, int order) 
+{
+   Mesh *mesh = fes.GetMesh();
+   IntegrationRule * ir;
+   int ir_order = 2*order+2;
+   IntegrationRules IntRule(0, Quadrature1D::GaussLobatto);
+   // IntegrationRules IntRule(0, Quadrature1D::GaussLegendre);
+
+   int dim = mesh->Dimension();
+   switch (dim)
+   {
+   case 1: 
+      ir = new IntegrationRule(IntRules.Get(mfem::Geometry::SEGMENT, ir_order));
+      break;
+   case 2: 
+      ir = new IntegrationRule(IntRules.Get(mfem::Geometry::SQUARE, ir_order));
+      break;   
+   default:
+      ir = new IntegrationRule(IntRules.Get(mfem::Geometry::CUBE, ir_order));
+      break;
+   }
+   return ir;
+}
+IntegrationRule * TensorIntegrationRule(int dim, int order) 
+{
+   IntegrationRule * ir;
+   int ir_order = 2*order+2;
+   // IntegrationRules IntRule(0, Quadrature1D::GaussLegendre);
+   IntegrationRules IntRule(0, Quadrature1D::GaussLobatto);
+   switch (dim)
+   {
+   case 1: 
+      ir = new IntegrationRule(IntRules.Get(mfem::Geometry::SEGMENT, ir_order));
+      break;
+   case 2: 
+      ir = new IntegrationRule(IntRules.Get(mfem::Geometry::SQUARE, ir_order));
+      break;   
+   default:
+      ir = new IntegrationRule(IntRules.Get(mfem::Geometry::CUBE, ir_order));
+      break;
+   }
+   return ir;
+}
+
+void KronMult(const Vector & x, const Vector & y, Vector & z)
+{
+   int n = x.Size();
+   int m = y.Size();
+   z.SetSize(n*m);
+   for (int i=0; i<n; i++)
+   {
+      for (int j = 0; j<m; j++)
+      {
+         z(i*m+j) = x(i)*y(j);
+      }
+   }
+}
+
+void KronMult(const Vector & x, const Vector & y, const Vector & z, Vector & w)
+{
+   Vector xy;
+   KronMult(x,y,xy);
+   KronMult(xy,z,w);
+}
+
+void AlterLS(DenseMatrix & T, Vector & vecA, Vector & vecB)
+{
+   int n = T.Height();
+   int m = T.Width();
+   Vector x(n);
+   Vector temp(m);
+   double s;
+   Vector y(m); y.Randomize(1); y /= y.Norml2();
+   int maxit = 3;
+   for (int i=0; i<maxit; i++)
+   {
+      T.Mult(y,x); x /= x.Norml2();
+      T.MultTranspose(x,temp); y=temp; y/= y.Norml2();
+      s = InnerProduct(y,temp);
+   }
+   vecA = x; vecA *= sqrt(s);
+   vecB = y; vecB *= sqrt(s);
+}
+
+void AlterLS(DenseTensor & T, Vector & vecA, Vector & vecB, Vector & vecC)
+{
+   int n = T.SizeI();
+   int m = T.SizeJ();
+   int l = T.SizeK();
+
+   DenseMatrix A0(n,l*m);
+   DenseMatrix A1(m,l*n);
+   DenseMatrix A2(l,n*m);
+
+   for (int i=0; i<n; ++i)
+   {
+      for (int j=0; j<m; ++j)
+      {
+         for (int k=0; k<l; ++k)
+         {
+            A0(i,j+k*m) = T(i,j,k);
+            A1(j,i+k*n) = T(i,j,k);
+            A2(k,i+j*n) = T(i,j,k);
+         }
+      }
+   }
+
+   double s;
+   Vector x(n);
+   Vector temp;
+   // Vector temp;
+   Vector y(m); y.Randomize(1); y/=y.Norml2();
+   Vector z(l); z.Randomize(2); z/=z.Norml2();
+
+   int maxit = 3;
+   for (int i = 0; i<maxit; i++)
+   {
+      KronMult(z,y,temp); A0.Mult(temp,x); s=x.Norml2(); x/=s;
+      KronMult(z,x,temp); A1.Mult(temp,y); s=y.Norml2(); y/=s;
+      KronMult(y,x,temp); A2.Mult(temp,z); s=z.Norml2(); z/=s;
+   }
+   vecA = x; vecA *= cbrt(s);
+   vecB = y; vecB *= cbrt(s);
+   vecC = z; vecC *= cbrt(s);
+}
+
+ElementTPFunctionCoefficient::ElementTPFunctionCoefficient(FiniteElementSpace &fes, int iel, Coefficient &cf)
+{
+   coeff_avg = 0.0;
+   dim = fes.GetMesh()->Dimension();
+   ElementTransformation * Tr;
+   const IntegrationRule * ir;
+   const FiniteElement * fe = fes.GetFE(iel);
+   ir = TensorIntegrationRule(fes,fe->GetOrder()); 
+   int nint1D = (dim == 2) ? sqrt(ir->GetNPoints()) : cbrt(ir->GetNPoints());
+   int nrintx = nint1D;
+   int nrinty = nint1D;
+   int nrintz = (dim ==2) ? 0 : nint1D;
+   if (dim == 2)
+   {
+      A.SetSize(nrintx,nrinty);
+      Tr = fes.GetElementTransformation(iel);
+      nint = ir->GetNPoints();
+      for (int i = 0; i < nint; i++)
+      {
+         int nint1D = sqrt(ir->GetNPoints());
+         int iy = i/nint1D;
+         int ix = i%nint1D;
+         const IntegrationPoint &ip = ir->IntPoint(i);
+         Vector x(dim);
+         Tr->Transform(ip,x);
+         double val = cf.Eval(*Tr, ip);
+         A(ix,iy) = val;
+         coeff_avg += val;
+      }
+      AlterLS(A,VecX,VecY);
+      coeff_avg /= (double)nint;
+   }
+   else
+   {
+      T.SetSize(nrintx,nrinty,nrintz);
+      Tr = fes.GetElementTransformation(iel);
+      for (int i = 0; i < ir->GetNPoints(); i++)
+      {
+         int nint1D = cbrt(ir->GetNPoints());
+         int iz = i/(nint1D*nint1D);
+         int iy = (i - iz*nint1D*nint1D)/nint1D;
+         int ix = (i - iz*nint1D*nint1D)%nint1D;
+         const IntegrationPoint &ip = ir->IntPoint(i);
+         Tr->SetIntPoint (&ip);
+         double val = cf.Eval(*Tr, ip);
+         T(ix,iy,iz) = val;
+      }
+      AlterLS(T,VecX,VecY,VecZ);
+   }
+   delete ir;
+}
+
+
+double ElementTPFunctionCoefficient::Eval(ElementTransformation &T, const IntegrationPoint &ip)
+{
+   double val = 0.0;
+   switch (coord)
+   {
+      case 0: 
+      {
+         if (orient == 1)
+         {
+            val = VecX(nintx++);
+         }
+         else
+         {
+            int nend = VecX.Size();   
+            val = VecX(nend-1-nintx++);
+         }
+      }
+      break;
+      case 1: 
+      {
+         if (orient == 1)
+         {
+            val = VecY(ninty++);
+         }
+         else
+         {
+            int nend = VecY.Size();   
+            val = VecY(nend-1-ninty++);
+         }
+      } 
+      break;
+      case 2: 
+      {
+         if (orient == 1)
+         {
+            val = VecZ(nintz++); 
+         }
+         else
+         {
+            int nend = VecZ.Size();   
+            val = VecZ(nend-1-nintz++);
+         }
+      }
+      break;
+      case -1:
+      {
+         int nint1D = sqrt(nint);
+         int iy = mint/nint1D;
+         int ix = mint%nint1D;
+         val = VecX(ix) * VecY(iy);
+         mint++;
+      }
+      break;
+      case -2:
+      {
+         int nint1D = cbrt(nint);
+         int iz = mint/(nint1D*nint1D);
+         int iy = (mint - iz*nint1D*nint1D)/nint1D;
+         int ix = (mint - iz*nint1D*nint1D)%nint1D;
+         val = VecX(ix) * VecY(iy) * VecZ(iz);
+         mint++;
+      }
+      break;
+      default: MFEM_ABORT("ElementTPFunctionCoefficient::Eval: Wrong coord choice"); 
+         break;
+   }
+   return val;
+}
+
+
+void TPElementTransformation::Setup2D()
+{
+   Mesh * mesh = fes->GetMesh();
+   MFEM_VERIFY(dim == 2, "Wrong dimension");
+   int nel = mesh->GetNE();
+   TransA1D.SetSize(nel,dim);
+   TransB1D.SetSize(nel,dim);
+   // Get ElementTransformations for the 2D elements
+   for (int iel = 0; iel <nel; iel++)
+   {
+      // allocate memory for Trans1D
+      for (int d = 0; d<dim; d++)
+      {
+         TransA1D[iel][d] = new Vector;
+         TransB1D[iel][d] = new Vector;
+      }
+      ElementTransformation * T = mesh->GetElementTransformation(iel);
+      // Populate integration points and get the K = adj(J)/sqrt(detJ);
+      // Store K_11^2 + K_12^2 
+      //       K_21^2 + K_22^2 
+
+      const FiniteElement * fe = fes->GetFE(iel);
+      const IntegrationRule * ir = TensorIntegrationRule(*fes,fe->GetOrder()); 
+      int nint = ir->GetNPoints();
+      int nint1D = sqrt(nint);
+      DenseMatrix A, B;
+      A.SetSize(nint1D,nint1D);
+      B.SetSize(nint1D,nint1D);
+
+      for (int i = 0; i < nint; i++)
+      {
+         int iy = i/nint1D;
+         int ix = i%nint1D;
+         const IntegrationPoint &ip = ir->IntPoint(i);
+         T->SetIntPoint(&ip);
+
+         double detJ = T->Weight();
+         const DenseMatrix & adjJ = T->AdjugateJacobian();
+         DenseMatrix JtJ(adjJ.Height());
+         MultAtB(adjJ,adjJ,JtJ);
+         JtJ *= 1.0/abs(detJ);
+         // DenseMatrix adjtt(adjJ.Height());
+         // MultAtB(adjJ, adjJ, adjtt);
+         double valA = pow(abs(adjJ(0,0))+abs(adjJ(0,1)),2)/abs(detJ);
+         // double valA = abs(JtJ(0,0))+abs(JtJ(0,1));
+         double valB = pow(abs(adjJ(1,1))+abs(adjJ(1,0)),2)/abs(detJ);
+         // double valB = abs(JtJ(1,1))+abs(JtJ(1,0));
+         A(ix,iy) = valA;
+         B(ix,iy) = valB;
+      }
+      AlterLS(A,*TransA1D[iel][0],*TransA1D[iel][1]);
+      AlterLS(B,*TransB1D[iel][0],*TransB1D[iel][1]);
+   }
+}
+void TPElementTransformation::Setup3D()
+{
+   Mesh * mesh = fes->GetMesh();
+   MFEM_VERIFY(dim == 3, "Wrong dimension");
+   int nel = mesh->GetNE();
+   TransA1D.SetSize(nel,dim);
+   TransB1D.SetSize(nel,dim);
+   TransC1D.SetSize(nel,dim);
+   // Get ElementTransformations for the 3D elements
+   for (int iel = 0; iel <nel; iel++)
+   {
+      // allocate memory for Trans1D
+      for (int d = 0; d<dim; d++)
+      {
+         TransA1D[iel][d] = new Vector;
+         TransB1D[iel][d] = new Vector;
+         TransC1D[iel][d] = new Vector;
+      }
+      ElementTransformation * T = mesh->GetElementTransformation(iel);
+      // Populate integrations points and get the K = adj(J)/sqrt(detJ);
+      // Store K_11^2 + K_12^2 + K_13^2
+      //       K_21^2 + K_22^2 + K_23^2
+      //       K_31^2 + K_32^2 + K_33^2
+
+      const FiniteElement * fe = fes->GetFE(iel);
+      const IntegrationRule * ir = TensorIntegrationRule(*fes,fe->GetOrder()); 
+      int nint = ir->GetNPoints();
+      int nint1D = cbrt(nint);
+      DenseTensor A, B, C;
+      A.SetSize(nint1D,nint1D,nint1D);
+      B.SetSize(nint1D,nint1D,nint1D);
+      C.SetSize(nint1D,nint1D,nint1D);
+
+      for (int i = 0; i < nint; i++)
+      {
+         int iz = i/(nint1D*nint1D);
+         int iy = (i - iz*nint1D*nint1D)/nint1D;
+         int ix = (i - iz*nint1D*nint1D)%nint1D;
+
+         const IntegrationPoint &ip = ir->IntPoint(i);
+         T->SetIntPoint(&ip);
+
+         double detJ = T->Weight();
+         const DenseMatrix & adjJ = T->AdjugateJacobian();
+         // DenseMatrix JtJ(adjJ.Height());
+         // MultAtB(adjJ,adjJ,JtJ);
+         // JtJ *= 1.0/abs(detJ);
+         // double valA = abs(JtJ(0,0))+abs(JtJ(0,1))+abs(JtJ(0,2));
+         double valA = pow(abs(adjJ(0,0))+abs(adjJ(0,1))+abs(adjJ(0,2)),2)/abs(detJ);
+         // double valB = abs(JtJ(1,0))+abs(JtJ(1,1))+abs(JtJ(1,2));
+         double valB = pow(abs(adjJ(1,0))+abs(adjJ(1,1))+abs(adjJ(1,2)),2)/abs(detJ);
+         // double valC = abs(JtJ(2,0))+abs(JtJ(2,1))+abs(JtJ(2,2));
+         double valC = pow(abs(adjJ(2,0))+abs(adjJ(2,1))+abs(adjJ(2,2)),2)/abs(detJ);
+         A(ix,iy,iz) = valA;
+         B(ix,iy,iz) = valB;
+         C(ix,iy,iz) = valC;
+      }
+      AlterLS(A,*TransA1D[iel][0],*TransA1D[iel][1],*TransA1D[iel][2]);
+      AlterLS(B,*TransB1D[iel][0],*TransB1D[iel][1],*TransB1D[iel][2]);
+      AlterLS(C,*TransC1D[iel][0],*TransC1D[iel][1],*TransC1D[iel][2]);
+      // for (int i = 0; i<dim; i++)
+      // {
+      //    cout << "TransA["<<i<<"] = "; TransA1D[iel][i]->Print(cout, TransA1D[iel][i]->Size()); 
+      //    cout << "TransB["<<i<<"] = "; TransB1D[iel][i]->Print(cout, TransB1D[iel][i]->Size()); 
+      //    cout << "TransC["<<i<<"] = "; TransC1D[iel][i]->Print(cout, TransC1D[iel][i]->Size()); 
+      // }
+      // cin.get();
+   }
+}
+
+TPElementTransformation::TPElementTransformation(FiniteElementSpace &fes_) 
+: fes(&fes_)
+{
+   dim = fes->GetMesh()->Dimension();
+   if (dim == 2) 
+   {
+      Setup2D();
+   }
+   else
+   {
+      Setup3D();
+   }
+} 
+
+
+void GetVertexToEdgeCount(const Mesh * mesh, DenseMatrix & edge_counts)
+{
+   // serial
+   int dim = mesh->Dimension();
+   int nv = mesh->GetNV();
+   int ne = mesh->GetNEdges();
+   int nel = mesh->GetNE();
+
+   Array<int> ibeg(dim), iend(dim), inc(dim);
+   // loop through axis
+   for (int axis = 0; axis<dim; axis++)
+   {
+      switch (axis)
+      {
+         case 0:  ibeg[axis] = 0; inc[axis] = 2; iend[axis] = (dim == 2) ? 4 : 8  ; break; // "x" edges
+         case 1:  ibeg[axis] = 1; inc[axis] = 2; iend[axis] = (dim == 2) ? 4 : 8  ; break; // "y" edges
+         case 2:  ibeg[axis] = 8; inc[axis] = 1; iend[axis] = 12; break;                   // "z edges"
+         default: MFEM_ABORT("This should be unreachable"); break; 
+      }
+   }
+
+   Array<bool> edge_marker(ne); 
+   edge_counts.SetSize(nv,dim);
+   edge_counts = 0.0;
+   Array<int> edge_owned;
+   bool par = false;
+#ifdef MFEM_USE_MPI
+   const ParMesh * pmesh = dynamic_cast<const ParMesh *>(mesh);
+   if (pmesh) par = true;
+#endif   
+   edge_owned.SetSize(ne); edge_owned = 0;
+   ND_FECollection fec(1, dim);
+   if (par) 
+   {
+      ParFiniteElementSpace aux_fes(const_cast<ParMesh *>(pmesh), &fec);
+      int mytoffset = aux_fes.GetMyTDofOffset();
+      int tsize = aux_fes.GetTrueVSize();
+      Array<int> dofs;
+      for (int i=0; i<ne; ++i)
+      {
+         aux_fes.GetEdgeDofs(i, dofs);
+         const int ldof = (dofs[0] >= 0) ? dofs[0] : -1 - dofs[0];
+         int sign = aux_fes.GetLocalTDofNumber(ldof);
+         if (sign == -1) continue; // just a hack for now to work with AMR, need to rethink this.
+         int gdof = aux_fes.GetGlobalTDofNumber(ldof);
+         if (gdof >= mytoffset && gdof < mytoffset+tsize)
+         {
+            edge_owned[i] = 1;
+         }
+      }
+   }
+   else
+   {
+      edge_owned = 1;
+      if (mesh->Nonconforming())
+      {
+         NCMesh * ncmesh = mesh->ncmesh;
+         const mfem::NCMesh::NCList &nclist = ncmesh->GetEdgeList();
+         Array<mfem::NCMesh::Slave> ncslaves = nclist.slaves;
+         int ns = ncslaves.Size();
+         for (int i=0; i<ns; ++i)
+         {
+             int j = ncslaves[i].index;
+            edge_owned[j] = 0;
+         }
+      }
+   }
+
+   for (int d=0; d<dim; d++)
+   {
+      edge_marker = false;
+      for (int iel = 0; iel<nel; iel++)
+      {
+         Array<int> elem_edges, cor;
+         mesh->GetElementEdges(iel,elem_edges,cor);
+         for (int i = ibeg[d]; i<iend[d]; i += inc[d])
+         {
+            int edge = elem_edges[i];
+            if (edge_owned.Size())
+            {
+               // skip contributions from edge not owned by the proc
+               if (!edge_owned[edge]) continue; 
+            }
+            if (edge_marker[edge]) continue;
+            Array<int>edge_vert;
+            mesh->GetEdgeVertices(edge,edge_vert);
+            for (int j = 0; j < 2; j++)
+            {
+               int vertex = edge_vert[j];
+               edge_counts(vertex,d) += 1.0;
+            }
+            edge_marker[edge] = true;
+         }
+      }
+   }   
+}
+
+void GetDiffusionEdgeMatrix(int iedge, FiniteElementSpace * fes, 
+                               Vector & Jac1D, Vector & Coeff1D,
+                               const IntegrationRule *ir,  
+                               DenseMatrix &elmat, int orient)
+{
+   const FiniteElement * el = fes->GetEdgeElement(iedge);
+   int nd = el->GetDof();
+   int dim = el->GetDim();
+
+   DenseMatrix dshape(nd,dim);
+   elmat.SetSize(nd);
+   elmat = 0.0;
+   int nint = ir->GetNPoints();
+   
+   for (int i = 0; i < nint; i++)
+   {
+      const IntegrationPoint &ip = ir->IntPoint(i);
+      double w = ip.weight;
+      el->CalcDShape(ip, dshape);
+      int j = orient == -1 ? nint-i-1 : i;
+      double val = Coeff1D(j) * Jac1D(i);
+      w *= val;
+      AddMult_a_AAt(w, dshape, elmat);
+   }
+}
+
+void GetMassEdgeMatrix(int iedge, FiniteElementSpace * fes, 
+                       Vector & Jac1D, Vector & Coeff1D,
+                       const IntegrationRule *ir,  
+                       DenseMatrix &elmat, int orient)
+{
+   const FiniteElement * el = fes->GetEdgeElement(iedge);
+   int nd = el->GetDof();
+   Vector shape(nd);
+   elmat.SetSize(nd);
+   elmat = 0.0;
+   int nint = ir->GetNPoints();
+   for (int i = 0; i < nint; i++)
+   {
+      const IntegrationPoint &ip = ir->IntPoint(i);
+      double w = ip.weight;
+      el->CalcShape(ip, shape);
+      int j = orient == -1 ? nint-i-1 : i;
+      double val = Coeff1D(j) * Jac1D(i);
+      w *= val;
+      AddMult_a_VVt(w, shape, elmat);
+   }
+}
+
+void Get1DMatrices(FiniteElementSpace * fes, int iedge, int orient, 
+                               Vector & JacL, Vector & JacM, Vector & Coeff1D,
+                               const IntegrationRule *ir,  
+                               DenseMatrix &L, DenseMatrix & M)
+{
+   const FiniteElement * el = fes->GetEdgeElement(iedge);
+   int dim = el->GetDim();
+   int nd = el->GetDof();
+
+   DenseMatrix dshape(nd,dim);
+   Vector shape(nd);
+   L.SetSize(nd); L = 0.0;
+   M.SetSize(nd); M = 0.0;
+   int nint = ir->GetNPoints();
+   for (int i = 0; i < nint; i++)
+   {
+      const IntegrationPoint &ip = ir->IntPoint(i);
+      double w = ip.weight;
+      el->CalcDShape(ip, dshape);
+      el->CalcShape(ip, shape);
+      int j = orient == -1 ? nint-i-1 : i;
+      double wL = w*Coeff1D(j) * JacL(j);
+      double wM = w*Coeff1D(j) * JacM(j);
+      AddMult_a_AAt(wL, dshape, L);
+      AddMult_a_VVt(wM, shape, M);
+   } 
+}                               
+
+
+
+void TensorProductEssentialDofsMaps(const Array<int> & ess_tdof_list, 
+                                    const ParFiniteElementSpace * fes, 
+                                    Array<Array<int> *> & tmap, // local edge map
+                                    Array<Array<int>* > & non_ess_dofs) // element map
+{
+
+   MPI_Comm comm = fes->GetComm();
+   int num_procs,myid;
+   MPI_Comm_size(comm, &num_procs);
+   MPI_Comm_rank(comm, &myid);
+
+   // 1. Find the element local dofs that are essential and then identify 
+   //    edge vertices (in local numbering) that are produced from
+   //    Gather the vertex dofs to be eliminated for each edge in tmap.
+   // 2. Reconstruct the essential dofs for each element (this list might be different 
+   //    for each element e.g. l-shape/fichera mesh, where to keep the kronecker product 
+   //    structure an essential dof is not eliminated, see vertex (v) below)
+   //   .  .  .
+   //   .  .  .
+   //   .  .  v  .  .
+   //   .  .  .  .  .
+   //   .  .  .  .  .
+   // ----------------------------------------------------------------------
+   //
+   ParMesh * pmesh = fes->GetParMesh();
+   int dim = pmesh->Dimension();
+   int nredges = pmesh->GetNEdges();
+   int tsize = fes->GetTrueVSize();
+   int vsize = fes->GetVSize();
+   Vector tess_tdof_marker(tsize); tess_tdof_marker = 0.0;
+   Vector ess_tdof_marker(vsize); 
+   
+   for (int i = 0; i<ess_tdof_list.Size(); i++)
+   {
+      int tdof = ess_tdof_list[i];
+
+      tess_tdof_marker[tdof] = 1.0;
+   }
+
+   fes->GetProlongationMatrix()->Mult(tess_tdof_marker, ess_tdof_marker);
+
+   tmap.SetSize(nredges);
+   for (int i = 0; i<nredges; i++) { tmap[i] = new Array<int>(0); } 
+   int nel = pmesh->GetNE();
+   non_ess_dofs.SetSize(nel);
+   for (int iel = 0; iel<nel; iel++)
+   {
+      non_ess_dofs[iel] = new Array<int>(0);
+      Array<int> local_dofs; 
+      Array<int> local_tdofs; 
+      const FiniteElement &fe = *fes->GetFE(iel);
+      // mfem to Tensor basis map
+      const Array<int> &dmap =
+      dynamic_cast<const TensorBasisElement&>(fe).GetDofMap();
+      Array<int> dmapt(dmap.Size());
+      for (int i = 0; i<dmapt.Size(); i++)
+      {
+         dmapt[dmap[i]] = i;
+      }
+      Array<int> elem_dofs;
+      fes->GetElementDofs(iel,elem_dofs);
+      // get local index of ess_dofs
+      int n = elem_dofs.Size();
+      // loop through the vertices
+      for (int i = 0; i< elem_dofs.Size(); i++)
+      {
+         int ldof = elem_dofs[i];
+         if (!ess_tdof_marker[ldof]) // if not essential dof
+         {
+            local_dofs.Append(dmapt[i]); // append in local element dofs
+            local_tdofs.Append(ldof); // append in local element dofs
+         }
+      }
+      
+      // Find the possible vertex local dofs on the edges to be removed
+      int n1D = (dim == 2) ? sqrt(n) : cbrt(n);
+      Array<int> edges, cor;
+      pmesh->GetElementEdges(iel,edges,cor);
+      Array<int> eidx(dim);
+      eidx[0] = 0;
+      eidx[1] = 1;
+      if (dim == 3) eidx[2] = 8;
+      // mark edge local dofs
+      Array<Array<int> *> markers(dim);
+      for (int d = 0; d<dim; d++)
+      {
+         Array<int> marker(n1D); marker = 0;
+         for (int i = 0; i<local_dofs.Size(); i++)
+         {
+            int j = local_dofs[i];
+            int c = j/(n1D*n1D);
+            int l = c*n1D*n1D;
+            int k = (d == 0) ? (j-l)%n1D : d == 1 ? (j-l)/n1D : c; 
+            marker[k] = 1;
+         }
+         markers[d] = new Array<int>(marker);
+         // pick up the edge and orientation
+         int edge = edges[eidx[d]];
+         int orient = cor[eidx[d]];
+         const FiniteElement &fe = *fes->GetEdgeElement(edge);
+         // edge tensor product map
+         const Array<int> &emap =
+         dynamic_cast<const TensorBasisElement&>(fe).GetDofMap();
+         Array<int> edge_ldofs;
+         for (int i = 0; i<marker.Size(); i++)
+         {
+            if (!marker[i])
+            {
+               if (orient == 1)
+               {
+                  edge_ldofs.Append(emap[i]);
+               }
+               else
+               {
+                  int k = (emap[i] == 1) ? 0 : 1;
+                  edge_ldofs.Append(k);
+               }
+            }
+         }
+         edge_ldofs.Sort(); edge_ldofs.Unique();   
+         tmap[edge]->Append(edge_ldofs);
+         tmap[edge]->Sort();
+         tmap[edge]->Unique();
+      }
+      if (dim == 2)
+      {
+         for (int j = 0; j<n1D; j++)
+         {
+            if ((*markers[1])[j])
+            {
+               for (int i = 0; i<n1D; i++)
+               {
+                  if ((*markers[0])[i])
+                  {
+                     int ldof = n1D*j+i;
+                     non_ess_dofs[iel]->Append(ldof);
+                  }
+               }
+            }
+         }
+      }
+      else
+      {
+         for (int k = 0; k<n1D; k++)
+         {
+            if ((*markers[2])[k])
+            {
+               for (int j = 0; j<n1D; j++)
+               {
+                  if ((*markers[1])[j])
+                  {
+                     for (int i = 0; i<n1D; i++)
+                     {
+                        if ((*markers[0])[i])
+                        {
+                           int ldof = n1D*n1D*k + n1D*j+i;
+                           non_ess_dofs[iel]->Append(ldof);
+                        }
+                     }
+                  }
+               }
+            }
+         }
+      }
+      non_ess_dofs[iel]->Sort();
+      non_ess_dofs[iel]->Unique();
+      for (int i = 0; i<non_ess_dofs[iel]->Size(); i++)
+      {
+         int ldof = (*non_ess_dofs[iel])[i];
+         int tdof = elem_dofs[dmap[ldof]];
+         (*non_ess_dofs[iel])[i] = tdof;
+      }
+      for (int d=0; d<dim; d++)
+      {
+         delete markers[d];
+      }
+   }
+};

examples/fd/smoother-util.hpp

@@ -0,0 +1,112 @@
+
+#pragma once
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+IntegrationRule * TensorIntegrationRule(const FiniteElementSpace & fes, int order);
+IntegrationRule * TensorIntegrationRule(int dim, int order);
+
+void KronMult(const Vector & x, const Vector & y, Vector & z);
+void KronMult(const Vector & x, const Vector & y, const Vector & z, Vector & w);
+
+
+void AlterLS(DenseMatrix & T, Vector & vecA, Vector & vecB);
+
+void AlterLS(DenseTensor & T, Vector & vecA, Vector & vecB, Vector & vecC);
+
+class ElementTPFunctionCoefficient : public Coefficient//
+{
+private:
+   int dim;
+   DenseMatrix A;
+   DenseTensor T;
+   double coeff_avg;
+   Vector VecX;
+   Vector VecY;
+   Vector VecZ;
+   int orient;
+   int nint; // total num of integrations points
+   int mint=0; // counter for all the integrations points
+   int nintx = 0; // counter for the x integrations points
+   int ninty = 0; // counter for the y integration points
+   int nintz = 0; // counter for the z integration points
+   int coord = 0; // (indication flag for x,y or z coordinate)
+
+public:
+   ElementTPFunctionCoefficient(FiniteElementSpace &fes, int iel, Coefficient &cf);
+   double GetCoeffAvg() {return coeff_avg;}
+   void ResetCounters() { mint = nintx = ninty = nintz = 0; }
+   void ResetCounter(int c) 
+   { 
+      switch (c)
+      {
+         case 0: nintx = 0; break;
+         case 1: ninty = 0; break;
+         case 2: nintz = 0; break;
+         default: mint = 0; break;
+      }
+   }
+   void SetCoord(int coord_) { coord = coord_; }
+   void SetOrient(int orient_) { orient = orient_; }
+   virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip);
+   Vector * GetVecX(){return &VecX;}
+   Vector * GetVecY(){return &VecY;}
+   Vector * GetVecZ(){return &VecZ;}
+   virtual ~ElementTPFunctionCoefficient() { }
+};
+
+
+class TPElementTransformation
+{
+private:
+   int dim;
+   FiniteElementSpace * fes = nullptr;
+   Array2D<Vector *> TransA1D;
+   Array2D<Vector *> TransB1D;
+   Array2D<Vector *> TransC1D;
+   void Setup2D();
+   void Setup3D();
+public:
+   TPElementTransformation(FiniteElementSpace &fes_);
+   Vector * GetTPTransformation(int iel, int coord, int which_coeff)
+   {
+      switch(which_coeff)
+      {
+         case 0: return TransA1D[iel][coord]; break;
+         case 1: return TransB1D[iel][coord]; break;
+         case 2: 
+         {
+            MFEM_VERIFY(dim == 3, "Wrong coeff for this dimension");
+            return TransC1D[iel][coord]; 
+            break;
+         }
+         default: MFEM_ABORT("Wrong coeff selection"); return 0; break; 
+      }
+   }
+   ~TPElementTransformation() { }
+};
+
+void GetVertexToEdgeCount(const Mesh * mesh, DenseMatrix & edge_counts);
+
+void Get1DMatrices(FiniteElementSpace * fes, int iedge, int orient, 
+                               Vector & JacL, Vector & JacM, Vector & Coeff1D,
+                               const IntegrationRule *ir,  
+                               DenseMatrix &L, DenseMatrix & M);
+
+void GetDiffusionEdgeMatrix(int iedge, FiniteElementSpace * fes, 
+                               Vector & Jac1D, Vector & Coeff1D,
+                               const IntegrationRule *ir,  
+                               DenseMatrix &elmat, int orient);
+void GetMassEdgeMatrix(int iedge, FiniteElementSpace * fes, 
+                       Vector & Jac1D, Vector & Coeff1D,
+                       const IntegrationRule *ir,  
+                       DenseMatrix &elmat, int orient);
+
+void TensorProductEssentialDofsMaps(const Array<int> & ess_tdof_list, 
+                                    const ParFiniteElementSpace * fes, 
+                                    Array<Array<int> *> & tmap, // local edge map
+                                    Array<Array<int>* > & non_ess_dofs); // element map

linalg/solvers.cpp

@@ -242,9 +242,9 @@ OperatorChebyshevSmoother::OperatorChebyshevSmoother(const Operator &oper_,
    max_eig_estimate(max_eig_estimate_),
    N(d.Size()),
    dinv(N),
-   diag(d),
+   diag(&d),
    coeffs(order),
-   ess_tdof_list(ess_tdofs),
+   ess_tdof_list(&ess_tdofs),
    residual(N),
    oper(&oper_) { Setup(); }
 
@@ -263,13 +263,13 @@ OperatorChebyshevSmoother::OperatorChebyshevSmoother(const Operator &oper_,
      order(order_),
      N(d.Size()),
      dinv(N),
-     diag(d),
+     diag(&d),
      coeffs(order),
-     ess_tdof_list(ess_tdofs),
+     ess_tdof_list(&ess_tdofs),
      residual(N),
      oper(&oper_)
 {
-   OperatorJacobiSmoother invDiagOperator(diag, ess_tdofs, 1.0);
+   OperatorJacobiSmoother invDiagOperator(*diag, ess_tdofs, 1.0);
    ProductOperator diagPrecond(&invDiagOperator, oper, false, false);
 
 #ifdef MFEM_USE_MPI
@@ -306,16 +306,52 @@ OperatorChebyshevSmoother::OperatorChebyshevSmoother(const Operator* oper_,
                                power_tolerance) { }
 #endif
 
+#ifdef MFEM_USE_MPI
+OperatorChebyshevSmoother::OperatorChebyshevSmoother(const Operator &oper_,
+                                                     const Solver &prec_,
+                                                     int order_, MPI_Comm comm,
+                                                     int power_iterations, double power_tolerance)
+#else
+OperatorChebyshevSmoother::OperatorChebyshevSmoother(const Operator &oper_,
+                                                     const Solver &prec_,
+                                                     int order_, int power_iterations, double power_tolerance)
+#endif
+   : Solver(oper_.Height()),
+     order(order_),
+     diag(nullptr),
+     N(oper_.Height()),
+     coeffs(order),
+     ess_tdof_list(nullptr),
+     residual(N),
+     oper(&oper_),
+     prec(&prec_)
+{
+   ProductOperator Precond(prec, oper, false, false);
+
+#ifdef MFEM_USE_MPI
+   PowerMethod powerMethod(comm);
+#else
+   PowerMethod powerMethod;
+#endif
+   Vector ev(oper->Width());
+   max_eig_estimate = powerMethod.EstimateLargestEigenvalue(Precond, ev,
+                                                            power_iterations, power_tolerance);
+   Setup();
+}
+
+
 void OperatorChebyshevSmoother::Setup()
 {
-   // Invert diagonal
    residual.UseDevice(true);
-   auto D = diag.Read();
-   auto X = dinv.Write();
-   MFEM_FORALL(i, N, X[i] = 1.0 / D[i]; );
-   auto I = ess_tdof_list.Read();
-   MFEM_FORALL(i, ess_tdof_list.Size(), X[I[i]] = 1.0; );
-
+   // Invert diagonal
+   if (diag)
+   {
+      auto D = diag->Read();
+      auto X = dinv.Write();
+      auto I = ess_tdof_list->Read();
+      MFEM_FORALL(i, N, X[i] = 1.0 / D[i]; );
+      MFEM_FORALL(i, ess_tdof_list->Size(), X[I[i]] = 1.0; );
+   }
    // Set up Chebyshev coefficients
    // For reference, see e.g., Parallel multigrid smoothing: polynomial versus
    // Gauss-Seidel by Adams et al.
@@ -410,11 +446,22 @@ void OperatorChebyshevSmoother::Mult(const Vector& x, Vector &y) const
          residual = helperVector;
       }
 
-      // Scale residual by inverse diagonal
+      // Scale residual by inverse diagonal or apply the given preconditioner
       const int n = N;
-      auto Dinv = dinv.Read();
+
       auto R = residual.ReadWrite();
-      MFEM_FORALL(i, n, R[i] *= Dinv[i]; );
+      if (prec)
+      {
+         // No device yet
+         Vector z(residual.Size()); z = 0.0;
+         prec->Mult(residual,z);
+         residual = z;
+      }
+      else
+      {
+         auto Dinv = dinv.Read();
+         MFEM_FORALL(i, n, R[i] *= Dinv[i]; );
+      }
 
       // Add weighted contribution to y
       auto Y = y.ReadWrite();

linalg/solvers.hpp

@@ -239,12 +239,25 @@ public:
                              int order, MPI_Comm comm = MPI_COMM_NULL,
                              int power_iterations = 10,
                              double power_tolerance = 1e-8);
+
+   /** Chebyshev accelaration for the given preconditioner @a prec.
+       The largest eigenvalue of the  preconditoned operator
+       is estimated internally via a power method. The
+       accuracy of the estimated eigenvalue may be controlled via
+       power_iterations and power_tolerance. */
+   OperatorChebyshevSmoother(const Operator &oper_, const Solver &prec,
+                             int order, MPI_Comm comm = MPI_COMM_NULL,
+                             int power_iterations = 10,
+                             double power_tolerance = 1e-8);
 #else
    OperatorChebyshevSmoother(const Operator &oper_, const Vector &d,
                              const Array<int>& ess_tdof_list,
                              int order, int power_iterations = 10,
                              double power_tolerance = 1e-8);
 
+   OperatorChebyshevSmoother(const Operator &oper_, const Solver &prec,
+                             int order, int power_iterations = 10,
+                             double power_tolerance = 1e-8);
    /// Deprecated: see pass-by-reference version above
    MFEM_DEPRECATED
    OperatorChebyshevSmoother(const Operator* oper_, const Vector &d,
@@ -271,12 +284,13 @@ private:
    double max_eig_estimate;
    const int N;
    Vector dinv;
-   const Vector &diag;
+   const Vector * diag = nullptr;
    Array<double> coeffs;
-   const Array<int>& ess_tdof_list;
+   const Array<int> * ess_tdof_list = nullptr;
    mutable Vector residual;
    mutable Vector helperVector;
    const Operator* oper;
+   const Solver* prec=nullptr;
 };