testing amr results for l-shape domain problem

minor
adding AMR example for element-smoother
2021-08-17 19:52:22 -07:00 · 2021-08-04 17:45:20 -07:00 · 2021-07-23 17:08:19 -07:00 · 2021-07-22 16:40:48 -07:00 · 2021-07-22 12:53:29 -07:00 · 2021-07-22 12:53:02 -07:00
22 changed files with 3270 additions and 29 deletions
@@ -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;
+}
@@ -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(){};
+};
@@ -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;
+}
+
+
+
@@ -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);
+      }
+   }
+}
+
+
@@ -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.;
+}
@@ -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);
@@ -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
@@ -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
@@ -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];
+      }
+   }
+};
@@ -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
@@ -17,6 +17,7 @@ list(APPEND SRCS
  complex_operator.cpp
  constraints.cpp
  densemat.cpp
+  fdsolver.cpp
  symmat.cpp
  handle.cpp
  matrix.cpp
@@ -39,6 +40,7 @@ list(APPEND HDRS
  dinvariants.hpp
  symmat.hpp
  dtensor.hpp
+  fdsolver.hpp
  handle.hpp
  invariants.hpp
  kernels.hpp
@@ -50,6 +50,10 @@ dsyevr_(char *JOBZ, char *RANGE, char *UPLO, int *N, double *A, int *LDA,
        double *W, double *Z, int *LDZ, int *ISUPPZ, double *WORK, int *LWORK,
        int *IWORK, int *LIWORK, int *INFO);
 extern "C" void
+dgeev_(const char * jobvl, const char * jobvr, int *n, double * A, int * lda,
+       double * wr, double * wl, double * vl, int * ldvl, double * vr, int * ldvr,
+       double * work, int * lwork, int * info);
+extern "C" void
 dsyev_(char *JOBZ, char *UPLO, int *N, double *A, int *LDA, double *W,
       double *WORK, int *LWORK, int *INFO);
 extern "C" void
@@ -2859,6 +2863,101 @@ void AddMult_a_VVt(const double a, const Vector &v, DenseMatrix &VVt)
   }
 }

+void KronProd(const DenseMatrix & A, const DenseMatrix & B, DenseMatrix & C)
+{
+   const int ah = A.Height();
+   const int aw = A.Width();
+   const int bh = B.Height();
+   const int bw = B.Width();
+
+   C.SetSize(ah*bh,aw*bw);
+   const double * ad = A.Data();
+   const double * bd = B.Data();
+   double * cd = C.Data();
+
+   for (int ja = 0; ja<aw; ++ja)
+      for (int jb = 0; jb<bw; ++jb)
+         for (int ia = 0; ia<ah; ++ia)
+            for (int ib = 0; ib<bh; ++ib)
+               cd[bh*ia + ib + ah*bh*(bw*ja + jb)]
+                  = ad[ia + ja * ah] * bd[ib + jb*bh];
+}
+
+void KronMult(const DenseMatrix &A, const DenseMatrix &B, const Vector &r,
+              Vector & z)
+{
+   const int nA = A.Height();
+   const int mA = A.Width();
+   const int nB = B.Height();
+   const int mB = B.Width();
+   const int nr = r.Size();
+   MFEM_VERIFY(nr == mA*mB, "Wrong size of Vector r");
+   z.SetSize(nA*nB);
+   DenseMatrix R(r.GetData(),mB,mA);
+   DenseMatrix X(nB,mA);
+   DenseMatrix Y(z.GetData(),nB,nA);
+   Mult(B,R,X);
+   MultABt(X,A,Y);
+}
+
+void KronMult(const DenseMatrix &A, const DenseMatrix &B, const DenseMatrix &R,
+              DenseMatrix & Z)
+{
+   const int nA = A.Height();
+   const int nB = B.Height();
+   const int nR = R.Height();
+   const int mR = R.Width();
+   Z.SetSize(nA*nB,mR);
+
+   Vector r,z;
+   double * dataR = R.Data();
+   for (int i = 0; i<mR; i++)
+   {
+      r.SetDataAndSize(&dataR[i*nR],nR);
+      KronMult(A,B,r,z);
+      Z.SetCol(i,z);
+   }
+}
+
+void KronMult(const DenseMatrix &A, const DenseMatrix &B, const DenseMatrix &C,
+              const Vector &r, Vector &z)
+{
+   const int nA = A.Height();
+   const int mA = A.Width();
+   const int nB = B.Height();
+   const int mB = B.Width();
+   const int nC = C.Height();
+   const int mC = C.Width();
+   const int nr = r.Size();
+   MFEM_VERIFY(nr == mA*mB*mC, "Wrong size of Vector r");
+   z.SetSize(nA*nB*nC);
+
+   double * dataR = r.GetData();
+   DenseMatrix R(dataR,mC,mA*mB);
+   DenseMatrix X(nC,mA*mB);
+   Mult(C,R,X);
+   X.Transpose();
+   DenseMatrix Z(z.GetData(),mA*mB,nC);
+   KronMult(A,B,X,Z);
+   Z.Transpose();
+}
+
+void KronMult(const Array<DenseMatrix *> & A, const Vector & r, Vector & z)
+{
+   int dim = A.Size();
+   if (dim == 2)
+   {
+      KronMult(*A[0],*A[1],r,z);
+   }
+   else if (dim == 3)
+   {
+      KronMult(*A[0],*A[1],*A[2], r,z);
+   }
+   else
+   {
+      MFEM_ABORT("KronMult::Wrong dimension");
+   }
+}

 bool LUFactors::Factor(int m, double TOL)
 {
@@ -3310,23 +3409,105 @@ DenseMatrixInverse::~DenseMatrixInverse()
   delete [] lu.ipiv;
 }

+void KronMult(const DenseMatrixInverse &A, const DenseMatrixInverse &B,
+              const Vector &r, Vector & z)
+{
+   // A and B are square matrices
+   z.SetSize(r.Size());
+   int nA = A.Height();
+   int nB = B.Height();
+   DenseMatrix R(r.GetData(),nB,nA);
+   DenseMatrix X(nB,nA);
+   B.Mult(R,X);
+   X.Transpose();
+   DenseMatrix Y(z.GetData(),nA,nB);
+   A.Mult(X,Y);
+   Y.Transpose();
+}

-DenseMatrixEigensystem::DenseMatrixEigensystem(DenseMatrix &m)
+void KronMult(const DenseMatrixInverse &A, const DenseMatrixInverse &B,
+              const DenseMatrix &R, DenseMatrix & Z)
+{
+   // A and B are square matrices
+   int nR = R.Height();
+   int mR = R.Width();
+   Z.SetSize(nR,mR);
+   Vector r(nR);
+   Vector z(nR);
+   double * dataR = R.GetData();
+   double * dataZ = Z.GetData();
+   for (int i = 0; i<mR; i++)
+   {
+      r.SetData(&dataR[i*nR]);
+      z.SetData(&dataZ[i*nR]);
+      KronMult(A,B,r,z);
+   }
+}
+
+void KronMult(const DenseMatrixInverse &A, const DenseMatrixInverse &B,
+              const DenseMatrixInverse &C, const Vector &r, Vector & z)
+{
+   // A, B and C are square matrices
+   int n = r.Size();
+   z.SetSize(n);
+   int nA = A.Height();
+   int nB = B.Height();
+   int nC = C.Height();
+   double * dataR = r.GetData();
+   DenseMatrix R(dataR,nC,nA*nB);
+   DenseMatrix X(nC,nA*nB);
+   C.Mult(R,X);
+   X.Transpose();
+   DenseMatrix Z(z.GetData(),nC,nA*nB);
+   KronMult(A,B,X,Z);
+   Z.Transpose();
+}
+
+void KronMult(const Array<DenseMatrixInverse *> & A, const Vector & r,
+              Vector & z)
+{
+   int dim = A.Size();
+   if (dim == 2)
+   {
+      KronMult(*A[0],*A[1],r,z);
+   }
+   else if (dim == 3)
+   {
+      KronMult(*A[0],*A[1],*A[2], r,z);
+   }
+   else
+   {
+      MFEM_ABORT("KronMult::Wrong dimension");
+   }
+}
+
+DenseMatrixEigensystem::DenseMatrixEigensystem(DenseMatrix &m, bool sym_)
   : mat(m)
 {
   n = mat.Width();
   EVal.SetSize(n);
+   EVali.SetSize(n);
   EVect.SetSize(n);
   ev.SetDataAndSize(NULL, n);

 #ifdef MFEM_USE_LAPACK
+   sym = sym_;
   jobz = 'V';
-   uplo = 'U';
   lwork = -1;
   double qwork;
-   dsyev_(&jobz, &uplo, &n, EVect.Data(), &n, EVal.GetData(),
-          &qwork, &lwork, &info);
-
+   if (sym)
+   {
+      uplo = 'U';
+      dsyev_(&jobz, &uplo, &n, EVect.Data(), &n, EVal.GetData(),
+             &qwork, &lwork, &info);
+   }
+   else
+   {
+      char jobvl = 'N';
+      int ldvl = 1;
+      dgeev_(&jobvl,&jobz,&n, mat.GetData(), &n, EVal.GetData(), EVali.GetData(),
+             nullptr, &ldvl, EVect.GetData(), &n, &qwork, &lwork, &info);
+   }
   lwork = (int) qwork;
   work = new double[lwork];
 #endif
@@ -3338,6 +3519,7 @@ DenseMatrixEigensystem::DenseMatrixEigensystem(
     n(other.n)
 {
 #ifdef MFEM_USE_LAPACK
+   sym = other.sym;
   jobz = other.jobz;
   uplo = other.uplo;
   lwork = other.lwork;
@@ -3356,13 +3538,24 @@ void DenseMatrixEigensystem::Eval()
 #endif

 #ifdef MFEM_USE_LAPACK
-   EVect = mat;
-   dsyev_(&jobz, &uplo, &n, EVect.Data(), &n, EVal.GetData(),
-          work, &lwork, &info);
-
+   if (sym)
+   {
+      EVect = mat;
+      dsyev_(&jobz, &uplo, &n, EVect.Data(), &n, EVal.GetData(),
+             work, &lwork, &info);
+   }
+   else
+   {
+      char jobvl = 'N';
+      int ldvl = 1;
+      DenseMatrix T = mat; // mat is overwritten by dgeev
+      dgeev_(&jobvl,&jobz,&n, T.GetData(), &n, EVal.GetData(), EVali.GetData(),
+             nullptr, &ldvl, EVect.GetData(), &n, work, &lwork, &info);
+   }
   if (info != 0)
   {
-      mfem::err << "DenseMatrixEigensystem::Eval(): DSYEV error code: "
+      string lpck = (sym) ? "DSYEV" : "DGEEV";
+      mfem::err << "DenseMatrixEigensystem::Eval(): " << lpck << "error code: "
                << info << endl;
      mfem_error();
   }
@@ -523,6 +523,22 @@ void AddMult_a_VWt(const double a, const Vector &v, const Vector &w,
 /// VVt += a * v v^t
 void AddMult_a_VVt(const double a, const Vector &v, DenseMatrix &VVt);

+/// C = A ⊗ B
+void KronProd(const DenseMatrix & A, const DenseMatrix & B, DenseMatrix & C);
+
+/// z = (A ⊗ B) r  = vec(B R A^T), where R := vec^-1 (r)
+void KronMult(const DenseMatrix &A, const DenseMatrix &B, const Vector &r,
+              Vector & z);
+
+/// z = (A ⊗ B) R
+void KronMult(const DenseMatrix &A, const DenseMatrix &B, const DenseMatrix &R,
+              DenseMatrix & Z);
+
+/// z = ( A ⊗ B ⊗ C ) r
+void KronMult(const DenseMatrix &A, const DenseMatrix &B, const DenseMatrix &C,
+              const Vector &r, Vector & z);
+
+void KronMult(const Array<DenseMatrix *> & A, const Vector & r, Vector & z);

 /** Class that can compute LU factorization of external data and perform various
    operations with the factored data. */
@@ -685,16 +701,33 @@ public:
   virtual ~DenseMatrixInverse();
 };

+/// z = (A^-1 ⊗ B^-1) r  = vec(B^-1 R A^-T), where R := vec^-1 (r)
+void KronMult(const DenseMatrixInverse &Ainv, const DenseMatrixInverse &Binv,
+              const Vector &r, Vector & z);
+
+/// z = (A^-1 ⊗ B^-1) R
+void KronMult(const DenseMatrixInverse &Ainv, const DenseMatrixInverse &Binv,
+              const DenseMatrix &R, DenseMatrix & Z);
+
+/// z = ( A^-1 ⊗ B^-1 ⊗ C^-1 ) r
+void KronMult(const DenseMatrixInverse &Ainv, const DenseMatrixInverse &Binv,
+              const DenseMatrixInverse &Cinv, const Vector &r, Vector & z);
+
+void KronMult(const Array<DenseMatrixInverse *> & A, const Vector & r,
+              Vector & z);

 class DenseMatrixEigensystem
 {
   DenseMatrix &mat;
   Vector      EVal;
+   // Possible non zero imaginary part of Eigenvalues
+   Vector      EVali;
   DenseMatrix EVect;
   Vector ev;
   int n;

 #ifdef MFEM_USE_LAPACK
+   bool sym;
   double *work;
   char jobz, uplo;
   int lwork, info;
@@ -702,10 +735,10 @@ class DenseMatrixEigensystem

 public:

-   DenseMatrixEigensystem(DenseMatrix &m);
+   DenseMatrixEigensystem(DenseMatrix &m, bool sym_ = false);
   DenseMatrixEigensystem(const DenseMatrixEigensystem &other);
   void Eval();
-   Vector &Eigenvalues() { return EVal; }
+   Vector &Eigenvalues(bool imag = false) { return imag ? EVali : EVal; }
   DenseMatrix &Eigenvectors() { return EVect; }
   double Eigenvalue(int i) { return EVal(i); }
   const Vector &Eigenvector(int i)
@@ -0,0 +1,143 @@
+// 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.
+
+#include "linalg.hpp"
+
+namespace mfem
+{
+
+void KronProdInvDiag(const Vector & a, const Vector & b, Vector & dinv)
+{
+   int n = a.Size(), m = b.Size();
+   dinv.SetSize(n*m);
+
+   for (int j = 0; j<m; j++)
+      for (int i = 0; i<n; i++)
+      {
+         dinv(i*m+j) = 1./(a(i) + b(j));
+      }
+}
+
+void KronProdInvDiag(const Vector & a, const Vector & b,
+                     const Vector & c, Vector & dinv)
+{
+   int n = a.Size(), m = b.Size(), l = c.Size();
+   dinv.SetSize(n*m*l);
+
+   for (int k = 0; k<l; k++)
+      for (int j = 0; j<m; j++)
+         for (int i = 0; i<n; i++)
+         {
+            dinv(i*m*l+j*l+k) = 1./(a(i) + b(j) + c(k));
+         }
+}
+
+void KronProdInvDiag(const Array<Vector *> & X, Vector & dinv)
+{
+   int dim = X.Size();
+   if (dim == 1)
+   {
+      int n = X[0]->Size();
+      dinv.SetSize(n);
+      for (int i = 0; i<n; i++) { dinv(i) = 1./(*X[0])(i); }
+   }
+   else if (dim == 2)
+   {
+      KronProdInvDiag(*X[0], *X[1], dinv);
+   }
+   else if (dim == 3)
+   {
+      KronProdInvDiag(*X[0], *X[1], *X[2], dinv);
+   }
+   else
+   {
+      MFEM_ABORT("KronProdInvDiag::Wrong dimension");
+   }
+}
+
+#ifdef MFEM_USE_LAPACK
+
+FDSolver::FDSolver(const Array<DenseMatrix *> & A,
+                   const Array<DenseMatrix *> & B)
+{
+   MFEM_ASSERT(A.Size() == B.Size(), "DenseFDSolver: Incompatible Dimensions");
+   dim = A.Size();
+
+   int solver_size = 1;
+   for (int i = 0; i<dim; i++)
+   {
+      MFEM_ASSERT(A[i]->Height() == A[i]->Width(),
+                  "DenseFDSolver: Matrix is not square");
+      MFEM_ASSERT(B[i]->Height() == B[i]->Width(),
+                  "DenseFDSolver: Matrix is not square");
+      MFEM_ASSERT(A[i]->Height() == B[i]->Height(),
+                  "DenseFDSolver: Matrices A and B have incompatible size");
+      solver_size *= A[i]->Height();
+   }
+   this->height = solver_size;
+   this->width  = solver_size;
+   if (solver_size) { Setup(A,B); }
+}
+
+void FDSolver::Setup(const Array<DenseMatrix *> & A,
+                     const Array<DenseMatrix *> & B)
+{
+   EigSystem.SetSize(dim);
+   eigv.SetSize(dim);
+   Array<Vector *> evalues(dim);
+   SQ.SetSize(dim);
+   DenseMatrix D;
+   for (int i = 0; i<dim; i++)
+   {
+      DenseMatrixInverse Minv(*B[i]);
+      Minv.Mult(*A[i],D);
+      EigSystem[i] = new DenseMatrixEigensystem(D);
+      EigSystem[i]->Eval();
+      evalues[i] = &EigSystem[i]->Eigenvalues();
+      eigv[i] = &EigSystem[i]->Eigenvectors();
+      DenseMatrixInverse Qinv(*eigv[i]);
+      DenseMatrix Sdinv;
+      Minv.GetInverseMatrix(Sdinv);
+      SQ[i] = new DenseMatrix;
+      Qinv.Mult(Sdinv,*SQ[i]);
+   }
+   KronProdInvDiag(evalues,dinv);
+}
+
+
+void FDSolver::Mult(const Vector & r,Vector & z) const
+{
+   MFEM_ASSERT(height == r.Size(),
+               "DenseFDSolver::Mult: Inconsistent vector size");
+   if (r.Size() == 0) { return; }
+   Vector rtemp;
+   KronMult(SQ,r,rtemp);
+   // 2. Diagonal solve;
+   rtemp *= dinv;
+   // 3. Modify RHS; z <-- (Q1 x Q2) rtemp
+   KronMult(eigv,rtemp,z);
+}
+
+FDSolver::~FDSolver()
+{
+   if (height)
+   {
+      for (int i=0; i<dim; i++)
+      {
+         delete SQ[i];
+         delete EigSystem[i];
+      }
+   }
+}
+
+#endif // MFEM_USE_LAPACK
+
+} // namespace mfem
@@ -0,0 +1,60 @@
+// 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.
+
+#ifndef MFEM_FDSOLVER
+#define MFEM_FDSOLVER
+
+#include "../config/config.hpp"
+#include "densemat.hpp"
+
+namespace mfem
+{
+
+/// Computes the inverse diagonal dinv = (a⊗I + I⊗b)^-1
+/// where a, b are diagonal matrices and I is the identity of the
+/// appropriate size
+void KronProdInvDiag(const Vector & a, const Vector & b, Vector & dinv);
+
+/// Computes the inverse diagonal dinv = (a⊗I⊗I + I⊗b⊗I + I⊗I⊗c)^-1
+/// where a, b, c are diagonal matrices and I is the identity of the
+/// appropriate size
+void KronProdInvDiag(const Vector & a, const Vector & b,
+                     const Vector & c, Vector & dinv);
+
+void KronProdInvDiag(const Array<Vector *> & X, Vector & dinv);
+
+#ifdef MFEM_USE_LAPACK
+
+/// In 2D it solves the system (A_0 ⊗ B_1 + B_0 ⊗ A_1) z = r
+/// In 3D it solves the system
+///     (A_0 ⊗ B_1 ⊗ B_2 + B_0 ⊗ A_1 ⊗ B_2 + B_0 ⊗ B_1 ⊗ A_2) z = r
+class FDSolver: public Solver
+{
+private:
+   int dim = 2;
+   Array<DenseMatrixEigensystem *> EigSystem;
+   Array<DenseMatrix *> eigv; // eigenvectors
+   Array<DenseMatrix *> SQ;
+   mutable Vector dinv;
+   void Setup(const Array<DenseMatrix *> & A, const Array<DenseMatrix *> & B);
+public:
+   FDSolver(const Array<DenseMatrix *> & A, const Array<DenseMatrix *> & B);
+   virtual void SetOperator(const Operator &op) {}
+   virtual void Mult(const Vector &r, Vector &z) const;
+   virtual ~FDSolver();
+};
+
+#endif // MFEM_USE_LAPACK
+
+} // mfem name space
+
+
+#endif // MFEM_FDSOLVER
@@ -31,6 +31,7 @@
 #include "invariants.hpp"
 #include "constraints.hpp"
 #include "auxiliary.hpp"
+#include "fdsolver.hpp"

 #ifdef MFEM_USE_AMGX
 #include "amgxsolver.hpp"
@@ -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();
@@ -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;
 };


@@ -28,6 +28,7 @@ set(UNIT_TESTS_SRCS
  linalg/test_complex_operator.cpp
  linalg/test_constrainedsolver.cpp
  linalg/test_direct_solvers.cpp
+  linalg/test_fdsolver.cpp
  linalg/test_hypre_ilu.cpp
  linalg/test_ilu.cpp
  linalg/test_matrix_block.cpp
@@ -0,0 +1,126 @@
+// 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.
+
+#include "mfem.hpp"
+#include "unit_tests.hpp"
+
+using namespace mfem;
+
+#ifdef MFEM_USE_LAPACK
+
+TEST_CASE("FDSolver",
+          "[FDSolver]")
+{
+   double tol = 1e-10;
+
+   // SPD matrices
+   DenseMatrix A0(
+   {
+      {
+         1.29919, 0.61256, 0.82545
+      },
+      {0.61256, 0.57891, 0.39662},
+      {0.82545, 0.39662, 0.57541}
+   });
+
+   DenseMatrix A1(
+   {
+      {0.748236, 0.701663, 0.607517, 0.236740},
+      {0.701663, 0.809316, 0.713186, 0.256070},
+      {0.607517, 0.713186, 0.794221, 0.233943},
+      {0.236740, 0.256070, 0.233943, 0.083129}
+   });
+
+   DenseMatrix B0(
+   {
+      {0.13483, 0.51389, 0.43052},
+      {0.51389, 2.26750, 1.86331},
+      {0.43052, 1.86331, 1.59869}
+   });
+
+   DenseMatrix B1(
+   {
+      {0.94177, 1.02400, 1.14743, 0.35723},
+      {1.02400, 1.79087, 1.78708, 0.78304},
+      {1.14743, 1.78708, 2.06259, 0.80837},
+      {0.35723, 0.78304, 0.80837, 1.01798}
+   });
+
+
+   SECTION("2D")
+   {
+      Array<DenseMatrix *> A(2), B(2);
+      A[0] = &A0;  A[1] = &A1;
+      B[0] = &B0;  B[1] = &B1;
+
+      Vector y(12); y.Randomize(1);
+      Vector x(12), diff(12);
+
+      FDSolver S(A,B);
+      S.Mult(y,x);
+
+      DenseMatrix C1, C;
+      KronProd(A0, B1, C1);
+      KronProd(B0, A1, C);
+
+      C.Add(1., C1);
+
+      DenseMatrixInverse Cinv(C);
+      Cinv.Mult(y,diff);
+
+      diff-=x;
+      REQUIRE(diff.Norml2() < tol);
+   }
+
+
+   SECTION("3D")
+   {
+      DenseMatrix A2(
+      {
+         {1.14593, 0.76119},
+         {0.76119, 0.78993}
+      });
+      DenseMatrix B2(
+      {
+         {0.88088, 0.37899},
+         {0.37899, 0.45096}
+      });
+      Array<DenseMatrix *> A(3), B(3);
+      A[0] = &A0;  A[1] = &A1; A[2] = &A2;
+      B[0] = &B0;  B[1] = &B1; B[2] = &B2;
+
+      Vector y(24); y.Randomize(1);
+      Vector x(24), diff(24);
+
+      FDSolver S(A,B);
+      S.Mult(y,x);
+
+      DenseMatrix Temp, C0, C1, C;
+
+      KronProd(A0, B1, Temp);
+      KronProd(Temp, B2, C0);
+      KronProd(B0, A1, Temp);
+      KronProd(Temp, B2, C1);
+      KronProd(B0, B1, Temp);
+      KronProd(Temp, A2, C);
+
+      C.Add(1.,C0);
+      C.Add(1.,C1);
+
+      DenseMatrixInverse Cinv(C);
+      Cinv.Mult(y,diff);
+
+      diff-=x;
+      REQUIRE(diff.Norml2() < tol);
+   }
+}
+
+#endif // if MFEM_USE_LAPACK
@@ -239,6 +239,188 @@ TEST_CASE("DenseMatrix A*B^T methods",
   }
 }

+TEST_CASE("KronMult methods",
+          "[DenseMatrix]")
+{
+   double tol = 1e-12;
+   int nA = 3, mA = 4;
+   int nB = 5, mB = 6;
+   DenseMatrix A(nA,mA);
+   DenseMatrix B(nB,mB);
+
+   for (int i = 0; i<nA; i++)
+      for (int j = 0; j<mA; j++)
+      {
+         A(i,j) = ((double)rand()/(double)RAND_MAX);
+      }
+
+   for (int i = 0; i<nB; i++)
+      for (int j = 0; j<mB; j++)
+      {
+         B(i,j) = ((double)rand()/(double)RAND_MAX);
+      }
+
+   DenseMatrix AB;
+   KronProd(A,B,AB);
+
+   // (A ⊗ B) r
+   SECTION("KronMultABr")
+   {
+      Vector r(mA*mB); r.Randomize();
+      MFEM_VERIFY(r.Size() == AB.Width(), "Check r size");
+      Vector z0(AB.Height());
+      AB.Mult(r,z0);
+
+      Vector z1;
+      KronMult(A,B,r,z1);
+      MFEM_VERIFY(z0.Size() == z1.Size(), "Check z1 size");
+      z0-=z1;
+      REQUIRE(z0.Norml2() < tol);
+   }
+   // (A ⊗ B) R
+   SECTION("KronMultABR")
+   {
+      int nR = mA*mB;
+      int mR = 7;
+      DenseMatrix R(nR, mR);
+      for (int i = 0; i<nR; i++)
+         for (int j = 0; j<mR; j++)
+         {
+            R(i,j) = ((double)rand()/(double)RAND_MAX);
+         }
+
+      DenseMatrix Z0(nA*nB,mR);
+      Mult(AB,R,Z0);
+
+      DenseMatrix Z1;
+      KronMult(A,B,R,Z1);
+      MFEM_VERIFY(Z0.Height() == Z1.Height() &&
+                  Z0.Width() == Z1.Width(), "Check z1 size");
+      Z0-=Z1;
+
+      REQUIRE(Z0.MaxMaxNorm() < tol);
+
+   }
+
+   // (A ⊗ B ⊗ C) r
+   SECTION("KronMultABCr")
+   {
+      int nC = 7, mC = 2;
+      DenseMatrix C(nC, mC);
+      for (int i = 0; i<nC; i++)
+         for (int j = 0; j<mC; j++)
+         {
+            C(i,j) = ((double)rand()/(double)RAND_MAX);
+         }
+
+      DenseMatrix ABC;
+      KronProd(AB,C,ABC);
+      Vector r(mA*mB*mC); r.Randomize();
+      MFEM_VERIFY(r.Size() == ABC.Width(), "Check r size");
+      Vector z0(nA*nB*nC);
+      ABC.Mult(r,z0);
+
+      Vector z1;
+      KronMult(A,B,C,r,z1);
+      MFEM_VERIFY(z0.Size() == z1.Size(), "Check z1 size");
+      z0-=z1;
+      REQUIRE(z0.Norml2() < tol);
+   }
+}
+
+TEST_CASE("KronMultInv methods",
+          "[DenseMatrixInverse]")
+{
+   double tol = 1e-12;
+   int nA = 3;
+   int nB = 2;
+   DenseMatrix A(
+   {
+      { 1.0,  0.2, 3.4},
+      {-2.0, -1.0, 3.1},
+      { 0.7,  1.4,-0.9}
+   });
+   DenseMatrix B(
+   {
+      {-10.1, 5.7},
+      {-3.0,  4.2}
+   });
+
+   DenseMatrixInverse Ainv(A);
+   DenseMatrixInverse Binv(B);
+
+   DenseMatrix AB;
+   KronProd(A,B,AB);
+
+   // (A^-1 ⊗ B^-1) r
+   SECTION("KronMultInvABr")
+   {
+
+      Vector r(nA*nB); r.Randomize();
+      MFEM_VERIFY(r.Size() == AB.Width(), "Check r size");
+      Vector z0(AB.Height());
+      DenseMatrixInverse ABinv(AB);
+      ABinv.Mult(r,z0);
+
+      Vector z1;
+      KronMult(Ainv,Binv,r,z1);
+      MFEM_VERIFY(z0.Size() == z1.Size(), "Check z1 size");
+      z0-=z1;
+      REQUIRE(z0.Norml2() < tol);
+   }
+
+   // (A^-1 ⊗ B^-1) R
+   SECTION("KronMultInvABR")
+   {
+      int nR = nA*nB;
+      int mR = 7;
+      DenseMatrix R(nR, mR);
+      for (int i = 0; i<nR; i++)
+         for (int j = 0; j<mR; j++)
+         {
+            R(i,j) = ((double)rand()/(double)RAND_MAX);
+         }
+
+      DenseMatrixInverse ABinv(AB);
+      DenseMatrix Z0(nA*nB,mR);
+      ABinv.Mult(R,Z0);
+
+      DenseMatrix Z1;
+      KronMult(Ainv,Binv,R,Z1);
+      MFEM_VERIFY(Z0.Height() == Z1.Height() &&
+                  Z0.Width() == Z1.Width(), "Check z1 size");
+      Z0-=Z1;
+
+      REQUIRE(Z0.MaxMaxNorm() < tol);
+   }
+
+   // (A^-1 ⊗ B^-1 ⊗ C^-1) r
+   SECTION("KronMultInvABCr")
+   {
+      int nC = 4;
+      DenseMatrix C(
+      {
+         {-2.1, 1.6, -3.4,  17.5},
+         {-7.1, 1.3, -7.5, -12.5},
+         { 0.5, 5.7, -6.0, -0.5},
+         { 9.2, 0.3, -1.4, -14.9}
+      });
+      DenseMatrix ABC;
+      KronProd(AB,C,ABC);
+      DenseMatrixInverse ABCInv(ABC);
+      Vector r(nA*nB*nC); r.Randomize();
+      MFEM_VERIFY(r.Size() == ABC.Width(), "Check r size");
+      Vector z0(nA*nB*nC);
+      ABCInv.Mult(r,z0);
+
+      DenseMatrixInverse Cinv(C);
+      Vector z1;
+      KronMult(Ainv,Binv,Cinv,r,z1);
+      MFEM_VERIFY(z0.Size() == z1.Size(), "Check z1 size");
+      z0-=z1;
+      REQUIRE(z0.Norml2() < tol);
+   }
+}

 TEST_CASE("LUFactors RightSolve", "[DenseMatrix]")
 {
@@ -314,6 +496,88 @@ TEST_CASE("DenseTensor LinearSolve methods",
   }
 }

+#ifdef MFEM_USE_LAPACK
+
+TEST_CASE("EigenSystem methods",
+          "[DenseMatrix]")
+{
+   double tol = 1e-12;
+   SECTION("SPD Matrix")
+   {
+      DenseMatrix A({{0.56806, 0.29211, 0.48315, 0.70024},
+         {0.29211, 0.85147, 0.68123, 0.70689},
+         {0.48315, 0.68123, 1.07229, 1.02681},
+         {0.70024, 0.70689, 1.02681, 1.15468}
+      });
+      DenseMatrix V, AV(4);
+      Vector Lambda;
+      for (bool sym: { false, true })
+      {
+         DenseMatrixEigensystem  EigA(A,sym);
+         EigA.Eval();
+         V = EigA.Eigenvectors();
+         Lambda = EigA.Eigenvalues();
+         Mult(A,V,AV);
+         V.RightScaling(Lambda);
+         AV -= V;
+         REQUIRE(AV.MaxMaxNorm() < tol);
+      }
+   }
+
+   SECTION("Indefinite Matrix")
+   {
+      DenseMatrix A({{0.486278, 0.041135, 0.480727, 0.616026},
+         {0.523599, 0.119827, 0.087808, 0.415241},
+         {0.214454, 0.661631, 0.909626, 0.744259},
+         {0.107007, 0.630604, 0.077862, 0.221006}
+      });
+      DenseMatrixEigensystem  EigA(A);
+      EigA.Eval();
+
+      Vector Lambda_r, Lambda_i;
+      // Real part of eigenvalues
+      Lambda_r = EigA.Eigenvalues();
+      // Imag part of eigenvalues
+      Lambda_i = EigA.Eigenvalues(true);
+
+      DenseMatrix V;
+      V = EigA.Eigenvectors();
+      // Real part of eigenvectors
+      DenseMatrix Vr(4), Vi(4);
+      Vr.SetCol(0,V.GetColumn(0));
+      Vr.SetCol(1,V.GetColumn(1));
+      Vr.SetCol(2,V.GetColumn(1));
+      Vr.SetCol(3,V.GetColumn(3));
+
+      // Imag part of eigenvectors
+      Vector vi(4); V.GetColumn(2,vi);
+      Vi.SetCol(0,0.);
+      Vi.SetCol(1,vi); vi *= -1.;
+      Vi.SetCol(2,vi);
+      Vi.SetCol(3,0.);
+
+      // Check that A*V = V * Lambda
+      // or A * (V_r + i V_i ) = (V_r + i V_i)*(Lamda_r + i Lambda_i)
+      // or  A * V_r = V_r * Lambda_r -  V_i * Lambda_i
+      // and A * V_i = V_r ( Lambda_i +  V_i * Lambda_r
+      DenseMatrix AVr(4), AVi(4);
+      Mult(A,Vr, AVr);
+      Mult(A,Vi, AVi);
+
+      DenseMatrix Vrlr = Vr; Vrlr.RightScaling(Lambda_r);
+      DenseMatrix Vrli = Vr; Vrli.RightScaling(Lambda_i);
+      DenseMatrix Vilr = Vi; Vilr.RightScaling(Lambda_r);
+      DenseMatrix Vili = Vi; Vili.RightScaling(Lambda_i);
+
+      AVr -= Vrlr; AVr+= Vili;
+      AVi -= Vrli; AVi-= Vilr;
+
+      REQUIRE(AVr.MaxMaxNorm() < tol);
+      REQUIRE(AVi.MaxMaxNorm() < tol);
+   }
+}
+
+#endif // if MFEM_USE_LAPACK
 TEST_CASE("DenseTensor copy", "[DenseMatrix][DenseTensor]")
 {
   DenseTensor t1(2,3,4);
Author	SHA1	Message	Date
psocratis	174b930544	testing amr results for l-shape domain problem	2021-08-17 19:52:22 -07:00
psocratis	afc10818a0	minor	2021-08-04 17:45:20 -07:00
psocratis	e9088d9b84	adding AMR example for element-smoother	2021-07-23 17:08:19 -07:00
psocratis	d62c094e97	adding diffusionMG class with element-smoother	2021-07-22 16:40:48 -07:00
psocratis	9541aa6026	diffusion par example testing the smoothers	2021-07-22 12:53:29 -07:00
psocratis	bf83b55578	adding new constructor to Chebychev smoother for a generic preconditioner (not just the diagonal)	2021-07-22 12:53:02 -07:00
psocratis	8a4d53098b	porting par-el-smoother 1st commit	2021-07-21 18:30:07 -07:00
psocratis	026d4834d2	Merge branch 'master' into fdsolver	2021-07-20 13:59:28 -07:00
psocratis	beeabc53c6	fixing failed test caused by auto merging with master	2021-06-11 10:56:12 -07:00
psocratis	f0124d2fc5	merge master	2021-06-11 10:48:55 -07:00
psocratis	f8a7405fb8	addressing some of the review comments	2021-05-13 20:23:36 -07:00
psocratis	76d384d6be	fixing unused variable warning	2021-04-28 20:22:52 -07:00
psocratis	76382289d0	Adding checks if MFEM_USE_LAPACK	2021-04-28 20:10:37 -07:00
psocratis	e0d52ab9f7	make style	2021-04-28 19:29:57 -07:00
psocratis	2214130993	fixing valgrind complaint in KronMult	2021-04-28 19:22:30 -07:00
psocratis	2b760e8a0a	Additional unit test for eigensystems of indefinite matrices	2021-04-28 18:48:02 -07:00
psocratis	ecfd6ff848	reordering the matrices in fdsolver	2021-04-28 18:46:12 -07:00
psocratis	f96929cc9e	Fixing issue where the given matrix given to the eigen system is overwritten	2021-04-28 18:42:57 -07:00
psocratis	dd4eb1a7b9	Corrected comment for the fdsolver description	2021-04-27 20:52:48 -07:00
psocratis	b2dce094f7	Adding unit tests for fdsolver	2021-04-27 20:50:27 -07:00
psocratis	e96797c7db	fixed small bug in lapack eig solver	2021-04-27 20:49:38 -07:00
psocratis	157f3e1a43	Cleaning up fdsolver	2021-04-27 19:08:57 -07:00
psocratis	c9c181d25b	adding entry-wise mult operator *=	2021-04-27 19:08:17 -07:00
psocratis	9b5b9ddcbb	adding KronMult for collection of DenseMatrices	2021-04-27 19:06:25 -07:00
psocratis	b328746f4a	make style	2021-04-26 20:24:54 -07:00
psocratis	da6f519346	Started fdsolver	2021-04-26 20:24:21 -07:00
psocratis	806595ccc2	Adding lapack EigenSystem for general dense matrides	2021-04-26 20:19:22 -07:00
psocratis	48183748ba	KronMult for DenseMatrixInverse and unit tests	2021-04-23 17:22:03 -07:00
psocratis	06ccc3cc29	adding KronMult and unit tests	2021-04-23 15:51:12 -07:00