Squashing commits on ho-trt-meshless

.appveyor.yml

@@ -26,9 +26,9 @@ install:
 - cd ..
 
 # Install hypre
-- ps: Start-FileDownload 'https://github.com/hypre-space/hypre/archive/V2-10-0b.tar.gz'
-- 7z x V2-10-0b.tar.gz -so | 7z x -si -ttar > nul
-- cd hypre-2-10-0b
+- ps: Start-FileDownload 'https://github.com/hypre-space/hypre/archive/v2.18.2.tar.gz'
+- 7z x v2.18.2.tar.gz -so | 7z x -si -ttar > nul
+- cd hypre-2.18.2\src
 - cmake -H. -Bbuild -DHYPRE_USING_FEI=OFF -DMPI_C_INCLUDE_PATH="C:\Program Files (x86)\Microsoft SDKs\MPI\Include" -DMPI_C_LIBRARIES="C:\Program Files (x86)\Microsoft SDKs\MPI\Lib\x86\msmpi.lib" -DMPI_CXX_LIBRARIES="C:\Program Files (x86)\Microsoft SDKs\MPI\Lib\x86\msmpi.lib" -DMPI_CXX_INCLUDE_PATH="C:\Program Files (x86)\Microsoft SDKs\MPI\Include"
 - cmake --build build
 - cmake --build build --target install

.travis.yml

@@ -291,18 +291,18 @@ install:
 
    # hypre
    - if [ $MPI == "YES" ]; then
-        if [ ! -e hypre-2.10.0b/src/hypre/lib/libHYPRE.a ]; then
-           wget https://computation.llnl.gov/project/linear_solvers/download/hypre-2.10.0b.tar.gz --no-check-certificate;
-           rm -rf hypre-2.10.0b;
-           tar xvzf hypre-2.10.0b.tar.gz;
-           cd hypre-2.10.0b/src;
+        if [ ! -e hypre-2.18.2/src/hypre/lib/libHYPRE.a ]; then
+           wget https://github.com/hypre-space/hypre/archive/v2.18.2.tar.gz --no-check-certificate;
+           rm -rf hypre-2.18.2;
+           tar xvzf v2.18.2.tar.gz;
+           cd hypre-2.18.2/src;
            ./configure --disable-fortran --without-fei CC=mpicc CXX=mpic++;
            make -j3;
            cd ../..;
         else
             echo "Reusing cached hypre-2.10.0b/";
         fi;
-        ln -s hypre-2.10.0b hypre;
+        ln -s hypre-2.18.2 hypre;
      else
         echo "Serial build, not using hypre";
      fi

examples/ex9p.cpp

@@ -58,6 +58,73 @@ double inflow_function(const Vector &x);
 // Mesh bounding box
 Vector bb_min, bb_max;
 
+struct AIR_parameters
+{
+   int blocksize;
+   int distanceR;
+   std::string prerelax;
+   std::string postrelax;
+   int interp_type;
+   int relax_type;
+   int coarsen_type;
+   double strength_tolC;
+   double strength_tolR;
+   double filter_tolR;
+   double filterA_tol;
+};
+
+class AIR_prec : public Solver
+{
+private:
+   const HypreParMatrix *A;
+   HypreParMatrix A_s;
+
+   // Preconditioner/solvers for A
+   HypreBoomerAMG *AIR_solver;
+   const AIR_parameters &AIR;
+   int blocksize;
+
+public:
+
+   AIR_prec(const AIR_parameters &_AIR) :
+      AIR_solver(NULL), AIR(_AIR)
+   {
+      blocksize = AIR.blocksize;
+   }
+
+   void SetOperator(const Operator &op)
+   {
+      A = dynamic_cast<const HypreParMatrix *>(&op);
+      delete AIR_solver;
+
+      // Scale A by block-diagonal inverse
+      BlockInvScal(A, &A_s, NULL, NULL, blocksize, 0);
+      AIR_solver = new HypreBoomerAMG(A_s);
+      AIR_solver->SetLAIROptions(AIR.distanceR, AIR.prerelax,
+                                 AIR.postrelax, AIR.strength_tolC,
+                                 AIR.strength_tolR, AIR.filter_tolR,
+                                 AIR.interp_type, AIR.relax_type,
+                                 AIR.filterA_tol, AIR.coarsen_type,
+                                 -1, 1);
+      AIR_solver->SetPrintLevel(0);
+      AIR_solver->SetMaxLevels(50);
+   }
+
+   virtual void Mult(const Vector &x, Vector &y) const
+   {
+      // scale the rhs by block inverse and solve system
+      HypreParVector z_s;
+      BlockInvScal(A, NULL, &x, &z_s, blocksize, 2);
+      AIR_solver->Mult(z_s, y);
+   }
+
+   ~AIR_prec()
+   {
+      BlockInvScal(NULL, NULL, NULL, NULL, 0, -1);
+      delete AIR_solver;
+   }
+};
+
 class DG_Solver : public Solver
 {
 private:
@@ -65,7 +132,8 @@ private:
    SparseMatrix M_diag;
    HypreParMatrix *A;
    GMRESSolver linear_solver;
-   BlockILU prec;
+   // BlockILU prec;
+   Solver *prec;
    double dt;
 public:
    DG_Solver(HypreParMatrix &M_, HypreParMatrix &K_, const FiniteElementSpace &fes)
@@ -73,20 +141,40 @@ public:
         K(K_),
         A(NULL),
         linear_solver(M.GetComm()),
-        prec(fes.GetFE(0)->GetDof(),
-             BlockILU::Reordering::MINIMUM_DISCARDED_FILL),
         dt(-1.0)
    {
+      prec = new BlockILU(fes.GetFE(0)->GetDof(),
+                          BlockILU::Reordering::MINIMUM_DISCARDED_FILL);
       linear_solver.iterative_mode = false;
       linear_solver.SetRelTol(1e-9);
       linear_solver.SetAbsTol(0.0);
       linear_solver.SetMaxIter(100);
       linear_solver.SetPrintLevel(0);
-      linear_solver.SetPreconditioner(prec);
+      linear_solver.SetPreconditioner(*prec);
 
       M.GetDiag(M_diag);
    }
 
+   DG_Solver(HypreParMatrix &M_, HypreParMatrix &K_, const FiniteElementSpace &fes,
+             const AIR_parameters &_AIR)
+      : M(M_),
+        K(K_),
+        A(NULL),
+        linear_solver(M.GetComm()),
+        dt(-1.0)
+   {
+      prec = new AIR_prec(_AIR);
+      linear_solver.iterative_mode = false;
+      linear_solver.SetRelTol(1e-9);
+      linear_solver.SetAbsTol(0.0);
+      linear_solver.SetMaxIter(100);
+      linear_solver.SetPrintLevel(0);
+      linear_solver.SetPreconditioner(*prec);
+
+      M.GetDiag(M_diag);
+   }
+
+
    void SetTimeStep(double dt_)
    {
       if (dt_ != dt)
@@ -115,10 +203,12 @@ public:
 
    ~DG_Solver()
    {
+      delete prec;
       delete A;
    }
 };
 
+
 /** A time-dependent operator for the right-hand side of the ODE. The DG weak
     form of du/dt = -v.grad(u) is M du/dt = K u + b, where M and K are the mass
     and advection matrices, and b describes the flow on the boundary. This can
@@ -137,6 +227,8 @@ private:
 
 public:
    FE_Evolution(ParBilinearForm &_M, ParBilinearForm &_K, const Vector &_b);
+   FE_Evolution(ParBilinearForm &_M, ParBilinearForm &_K, const Vector &_b,
+                const AIR_parameters &_AIR);
 
    virtual void Mult(const Vector &x, Vector &y) const;
    virtual void ImplicitSolve(const double dt, const Vector &x, Vector &k);
@@ -169,7 +261,11 @@ int main(int argc, char *argv[])
    bool paraview = false;
    bool binary = false;
    int vis_steps = 5;
-
+   int solver_type = 1;
+   int basis_type = BasisType::GaussLobatto;
+   AIR_parameters AIR0 = {-1, 1, "", "FA", 100, 10, 10,
+                          0.1, 0.01, 0.0, 1e-4
+                         };
    int precision = 8;
    cout.precision(precision);
 
@@ -199,6 +295,10 @@ int main(int argc, char *argv[])
                   "Final time; start time is 0.");
    args.AddOption(&dt, "-dt", "--time-step",
                   "Time step.");
+   args.AddOption(&basis_type, "-b", "--basis-type",
+                  "DG finite element basis type. 0 for Gauss-Leg., 1 for Gauss-Lob.");
+   args.AddOption(&solver_type, "-st", "--solver-type",
+                  "Solver for implicit solves. 0 for ILU, 1 for pAIR-AMG.");
    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                   "--no-visualization",
                   "Enable or disable GLVis visualization.");
@@ -291,7 +391,7 @@ int main(int argc, char *argv[])
 
    // 7. Define the parallel discontinuous DG finite element space on the
    //    parallel refined mesh of the given polynomial order.
-   DG_FECollection fec(order, dim, BasisType::GaussLobatto);
+   DG_FECollection fec(order, dim, basis_type);
    ParFiniteElementSpace *fes = new ParFiniteElementSpace(pmesh, &fec);
 
    HYPRE_Int global_vSize = fes->GlobalTrueVSize();
@@ -300,6 +400,9 @@ int main(int argc, char *argv[])
       cout << "Number of unknowns: " << global_vSize << endl;
    }
 
+   // Get DG blocksize for AIR
+   AIR0.blocksize = fes->GetFE(0)->GetDof();
+
    // 8. Set up and assemble the parallel bilinear and linear forms (and the
    //    parallel hypre matrices) corresponding to the DG discretization. The
    //    DGTraceIntegrator involves integrals over mesh interior faces.
@@ -427,11 +530,19 @@ int main(int argc, char *argv[])
    // 10. Define the time-dependent evolution operator describing the ODE
    //     right-hand side, and perform time-integration (looping over the time
    //     iterations, ti, with a time-step dt).
-   FE_Evolution adv(*m, *k, *B);
+   FE_Evolution *adv;
+   if (solver_type == 1)
+   {
+      adv = new FE_Evolution(*m, *k, *B, AIR0);
+   }
+   else
+   {
+      adv = new FE_Evolution(*m, *k, *B);
+   }
 
    double t = 0.0;
-   adv.SetTime(t);
-   ode_solver->Init(adv);
+   adv->SetTime(t);
+   ode_solver->Init(*adv);
 
    bool done = false;
    for (int ti = 0; !done; )
@@ -498,6 +609,7 @@ int main(int argc, char *argv[])
    delete ode_solver;
    delete pd;
    delete dc;
+   delete adv;
 
    MPI_Finalize();
    return 0;
@@ -551,6 +663,58 @@ FE_Evolution::FE_Evolution(ParBilinearForm &_M, ParBilinearForm &_K,
    M_solver.SetPrintLevel(0);
 }
 
+// Implementation of class FE_Evolution
+FE_Evolution::FE_Evolution(ParBilinearForm &_M, ParBilinearForm &_K,
+                           const Vector &_b, const AIR_parameters &_AIR)
+   : TimeDependentOperator(_M.Height()),
+     b(_b),
+     M_solver(_M.ParFESpace()->GetComm()),
+     z(_M.Height())
+{
+   bool pa = _M.GetAssemblyLevel()==AssemblyLevel::PARTIAL;
+
+   if (pa)
+   {
+      MFEM_ABORT("This combination of options is untested.");
+      M.Reset(&_M, false);
+      K.Reset(&_K, false);
+   }
+   else
+   {
+      M.Reset(_M.ParallelAssemble(), true);
+      K.Reset(_K.ParallelAssemble(), true);
+   }
+
+   M_solver.SetOperator(*M);
+
+   Array<int> ess_tdof_list;
+   if (pa)
+   {
+      M_prec = new OperatorJacobiSmoother(_M, ess_tdof_list);
+      dg_solver = NULL;
+   }
+   else
+   {
+      HypreParMatrix &M_mat = *M.As<HypreParMatrix>();
+      HypreParMatrix &K_mat = *K.As<HypreParMatrix>();
+      HypreSmoother *hypre_prec = new HypreSmoother(M_mat, HypreSmoother::Jacobi);
+      M_prec = hypre_prec;
+
+      dg_solver = new DG_Solver(M_mat, K_mat, *_M.FESpace(), _AIR);
+   }
+
+   M_solver.SetPreconditioner(*M_prec);
+   M_solver.iterative_mode = false;
+   M_solver.SetRelTol(1e-9);
+   M_solver.SetAbsTol(0.0);
+   M_solver.SetMaxIter(100);
+   M_solver.SetPrintLevel(0);
+}
+
+// Solve the equation:
+//    u_t = M^{-1}(Ku + b),
+// by solving associated linear system
+//    (M - dt*K) d = K*u + b
 void FE_Evolution::ImplicitSolve(const double dt, const Vector &x, Vector &k)
 {
    K->Mult(x, z);

fem/bilinearform.cpp

@@ -269,6 +269,32 @@ void BilinearForm::AddBdrFaceIntegrator(BilinearFormIntegrator *bfi,
    bfbfi_marker.Append(&bdr_marker);
 }
 
+// ADDED //
+void BilinearForm::ResetDomainIntegrator (BilinearFormIntegrator * bfi)
+{
+   dbfi.DeleteAll();
+   dbfi.Append (bfi);
+}
+
+void BilinearForm::ResetBoundaryIntegrator (BilinearFormIntegrator * bfi)
+{
+   bbfi.DeleteAll();
+   bbfi.Append (bfi);
+}
+
+void BilinearForm::ResetInteriorFaceIntegrator (BilinearFormIntegrator * bfi)
+{
+   fbfi.DeleteAll();
+   fbfi.Append (bfi);
+}
+
+void BilinearForm::ResetBdrFaceIntegrator (BilinearFormIntegrator * bfi)
+{
+   bfbfi.DeleteAll();
+   bfbfi.Append (bfi);
+}
+// ADDED //
+
 void BilinearForm::ComputeElementMatrix(int i, DenseMatrix &elmat)
 {
    if (element_matrices)

fem/bilinearform.hpp

@@ -311,6 +311,20 @@ public:
    void AddBdrFaceIntegrator(BilinearFormIntegrator *bfi,
                              Array<int> &bdr_marker);
 
+   // ADDED //
+   /// Resets to new Domain Integrator.
+   void ResetDomainIntegrator(BilinearFormIntegrator *bfi);
+
+   /// Resets to new Boundary Integrator.
+   void ResetBoundaryIntegrator(BilinearFormIntegrator *bfi);
+
+   /// Resets to new interior Face Integrator.
+   void ResetInteriorFaceIntegrator(BilinearFormIntegrator *bfi);
+
+   /// Resets to new boundary Face Integrator.
+   void ResetBdrFaceIntegrator(BilinearFormIntegrator *bfi);
+   // ADDED
+
    void operator=(const double a)
    {
       if (mat != NULL) { *mat = a; }

fem/bilininteg.cpp

@@ -2507,6 +2507,12 @@ void DGTraceIntegrator::AssembleFaceMatrix(const FiniteElement &el1,
 
    double un, a, b, w;
 
+// ADDED //
+#ifdef MFEM_THREAD_SAFE
+   Vector shape1, shape2;
+#endif
+// ADDED //
+
    dim = el1.GetDim();
    ndof1 = el1.GetDof();
    Vector vu(dim), nor(dim);

fem/eltrans.cpp

@@ -556,4 +556,62 @@ void IntegrationPointTransformation::Transform (const IntegrationRule &ir1,
    }
 }
 
+
+////////////////////////////// ADDED /////////////////////////////////
+int IsoparametricTransformation::TransformBack(const Vector &pt,
+                                               IntegrationPoint &ip,
+                                               IntegrationPoint &xip)
+{
+   const int    max_iter = 32;
+   const double  ref_tol = 1e-12;
+   const double phys_tol = 1e-12*pt.Normlinf();
+
+   const int dim = FElem->GetDim();
+   const int sdim = PointMat.Height();
+   const int geom = FElem->GetGeomType();
+  // IntegrationPoint xip, prev_xip;
+   IntegrationPoint prev_xip;
+   double xd[3], yd[3], dxd[3], Jid[9];
+   Vector x(xd, dim), y(yd, sdim), dx(dxd, dim);
+   DenseMatrix Jinv(Jid, dim, sdim);
+   bool hit_bdr = false, prev_hit_bdr;
+
+   // Use the center of the element as initial guess
+ //  xip = Geometries.GetCenter(geom);
+ //  xip.Get(xd, dim); // xip -> x
+
+   for (int it = 0; it < max_iter; it++)
+   {
+      // Newton iteration:    x := x + J(x)^{-1} [pt-F(x)]
+      // or when dim != sdim: x := x + [J^t.J]^{-1}.J^t [pt-F(x)]
+      Transform(xip, y);
+      subtract(pt, y, y); // y = pt-y
+      if (y.Normlinf() < phys_tol) { ip = xip; return 0; }
+      SetIntPoint(&xip);
+      CalcInverse(Jacobian(), Jinv);
+      Jinv.Mult(y, dx);
+      x += dx;
+      prev_xip = xip;
+      prev_hit_bdr = hit_bdr;
+      xip.Set(xd, dim); // x -> xip
+      // If xip is ouside project it on the boundary on the line segment
+      // between prev_xip and xip
+      hit_bdr = !Geometry::ProjectPoint(geom, prev_xip, xip);
+      if (dx.Normlinf() < ref_tol) { ip = xip; return 0; }
+      if (hit_bdr)
+      {
+         xip.Get(xd, dim); // xip -> x
+         if (prev_hit_bdr)
+         {
+            prev_xip.Get(dxd, dim); // prev_xip -> dx
+            subtract(x, dx, dx);    // dx = xip - prev_xip
+            if (dx.Normlinf() < ref_tol) { return 1; }
+         }
+      }
+   }
+   ip = xip;
+   return 2;
+}
+////////////////////////////// ADDED /////////////////////////////////
+
 }

fem/eltrans.hpp

@@ -107,6 +107,11 @@ public:
        transformations. */
    virtual int TransformBack(const Vector &pt, IntegrationPoint &ip) = 0;
 
+   // ADDED //
+   virtual int TransformBack(const Vector &, IntegrationPoint &, 
+                             IntegrationPoint &) = 0;
+   // ADDED //
+
    virtual ~ElementTransformation() { }
 };
 
@@ -338,6 +343,11 @@ public:
       return inv_tr.Transform(v, ip);
    }
 
+   // ADDED //
+   virtual int TransformBack(const Vector &pt, IntegrationPoint &ip,
+                             IntegrationPoint &xip);
+   // ADDED //
+
    virtual ~IsoparametricTransformation() { }
 };
 

fem/fe.hpp

@@ -214,7 +214,9 @@ public:
    {
       Pk, ///< Polynomials of order k
       Qk, ///< Tensor products of polynomials of order k
-      rQk ///< Refined tensor products of polynomials of order k
+      rQk, ///< Refined tensor products of polynomials of order k
+      RBF, ///< Radial basis functions
+      RK, ///< Reproducing kernels
    };
 };
 
@@ -1721,9 +1723,6 @@ private:
 
    static Array2D<int> binom;
 
-   static void CalcMono(const int p, const double x, double *u);
-   static void CalcMono(const int p, const double x, double *u, double *d);
-
    static void CalcChebyshev(const int p, const double x, double *u);
    static void CalcChebyshev(const int p, const double x, double *u, double *d);
    static void CalcChebyshev(const int p, const double x, double *u, double *d,
@@ -1812,6 +1811,9 @@ public:
    static void CalcBernstein(const int p, const double x, double *u, double *d)
    { CalcBinomTerms(p, x, 1. - x, u, d); }
 
+   static void CalcMono(const int p, const double x, double *u);
+   static void CalcMono(const int p, const double x, double *u, double *d);
+
    static void CalcLegendre(const int p, const double x, double *u);
    static void CalcLegendre(const int p, const double x, double *u, double *d);
 
@@ -3013,6 +3015,488 @@ public:
                              DenseMatrix &hessian) const;
 };
 
+class RBFFunction
+{
+public:
+   static const double GlobalRadius; // functions with r>=GR are considered global
+   RBFFunction() { };
+   virtual ~RBFFunction() { }
+
+   // The r is a normalized distance
+   virtual double BaseFunction(double r) const = 0;
+   virtual double BaseDerivative(double r) const = 0;
+   virtual double BaseDerivative2(double r) const = 0;
+
+   // The support radius, outside of which the function is zero
+   virtual double Radius() const { return GlobalRadius; }
+
+   // Does function have compact support?
+   virtual bool CompactSupport() const { return false; }
+   
+   // This makes the shape parameter consistent across kernels
+   virtual double HNorm() const = 0;
+};
+
+class GaussianRBF : public RBFFunction
+{
+   // hNorm minimizes integral of Gaussian minus Wendland kernel over r=0,1
+   static const double hNorm;
+public:
+   GaussianRBF() { };
+   virtual ~GaussianRBF() { }
+   
+   virtual double BaseFunction(double r) const;
+   virtual double BaseDerivative(double r) const;
+   virtual double BaseDerivative2(double r) const;
+   
+   virtual double HNorm() const { return hNorm; }
+};
+
+class MultiquadricRBF : public RBFFunction
+{
+   // Same as inverse multiquadric
+   static const double hNorm;
+public:
+   MultiquadricRBF() { };
+   virtual ~MultiquadricRBF() { }
+   
+   virtual double BaseFunction(double r) const;
+   virtual double BaseDerivative(double r) const;
+   virtual double BaseDerivative2(double r) const;
+   
+   virtual double HNorm() const { return hNorm; }
+};
+
+class InvMultiquadricRBF : public RBFFunction
+{
+   // hNorm minimizes integral of Gaussian minus InvMQ kernel over r=0,0.5
+   static const double hNorm;
+public:
+   InvMultiquadricRBF() { };
+   virtual ~InvMultiquadricRBF() { }
+   
+   virtual double BaseFunction(double r) const;
+   virtual double BaseDerivative(double r) const;
+   virtual double BaseDerivative2(double r) const;
+   
+   virtual double HNorm() const { return hNorm; }
+};
+
+// Shifted down to be exactly zero at radius
+class CompactGaussianRBF : public RBFFunction
+{
+   static const double hNorm;
+   const double radius;
+   double multK, shiftK;
+   
+public:
+   CompactGaussianRBF(const double rad = 5.0);
+   virtual ~CompactGaussianRBF() { }
+   
+   virtual double BaseFunction(double r) const;
+   virtual double BaseDerivative(double r) const;
+   virtual double BaseDerivative2(double r) const;
+
+   virtual double Radius() const { return radius; }
+   
+   virtual double HNorm() const { return hNorm; }
+   virtual bool CompactSupport() const { return true; }
+};
+
+// Truncated at radius
+class TruncatedGaussianRBF : public RBFFunction
+{
+   static const double hNorm;
+   const double radius;
+public:
+   TruncatedGaussianRBF(const double rad = 5.0)
+      : radius(rad) { }
+   virtual ~TruncatedGaussianRBF() { }
+   
+   virtual double BaseFunction(double r) const;
+   virtual double BaseDerivative(double r) const;
+   virtual double BaseDerivative2(double r) const;
+
+   virtual double Radius() const { return radius; }
+   
+   virtual double HNorm() const { return hNorm; }
+   virtual bool CompactSupport() const { return true; }
+};
+
+class Wendland11RBF : public RBFFunction
+{
+   static const double radius;
+   
+public:
+   Wendland11RBF() { }
+   virtual ~Wendland11RBF() { }
+   
+   virtual double BaseFunction(double r) const;
+   virtual double BaseDerivative(double r) const;
+   virtual double BaseDerivative2(double r) const;
+
+   virtual double Radius() const { return radius; }
+   
+   virtual double HNorm() const { return 1.0 / radius; }
+   virtual bool CompactSupport() const { return true; }
+};
+
+class Wendland31RBF : public RBFFunction
+{
+   static const double radius;
+   
+public:
+   Wendland31RBF() { };
+   virtual ~Wendland31RBF() { }
+   
+   virtual double BaseFunction(double r) const;
+   virtual double BaseDerivative(double r) const;
+   virtual double BaseDerivative2(double r) const;
+
+   virtual double Radius() const { return radius; }
+   
+   virtual double HNorm() const { return 1.0 / radius; }
+   virtual bool CompactSupport() const { return true; }
+};
+
+class Wendland33RBF : public RBFFunction
+{
+   static const double radius;
+   
+public:
+   Wendland33RBF() { };
+   virtual ~Wendland33RBF() { }
+   
+   virtual double BaseFunction(double r) const;
+   virtual double BaseDerivative(double r) const;
+   virtual double BaseDerivative2(double r) const;
+
+   virtual double Radius() const { return radius; }
+   
+   virtual double HNorm() const { return 1.0 / radius; }
+   virtual bool CompactSupport() const { return true; }
+};
+
+// Choose the type of RBF to use
+class RBFType
+{
+public:
+   enum
+   {
+      Gaussian = 0,
+      Multiquadric = 1,
+      InvMultiquadric = 2,
+      TruncatedGaussian = 3,
+      CompactGaussian = 4,
+      Wendland11 = 5,
+      Wendland31 = 6,
+      Wendland33 = 7,
+      NumRBFTypes = 8
+   };
+
+   // Return the requested RBF
+   static RBFFunction *GetRBF(const int rbfType)
+   {
+      switch (rbfType)
+      {
+      case RBFType::Gaussian:
+         return new GaussianRBF();
+      case RBFType::Multiquadric:
+         return new MultiquadricRBF();
+      case RBFType::InvMultiquadric:
+         return new InvMultiquadricRBF();
+      case RBFType::TruncatedGaussian:
+         return new TruncatedGaussianRBF();
+      case RBFType::CompactGaussian:
+         return new CompactGaussianRBF();
+      case RBFType::Wendland11:
+         return new Wendland11RBF();
+      case RBFType::Wendland31:
+         return new Wendland31RBF();
+      case RBFType::Wendland33:
+         return new Wendland33RBF();
+      }
+      MFEM_ABORT("unknown RBF type");
+      return NULL;
+   }
+   
+   // Abort if rbfType is invalid
+   static int Check(const int rbfType)
+   {
+      MFEM_VERIFY(0 <= rbfType && rbfType < NumRBFTypes,
+                  "unknown RBF type: " << rbfType);
+      return rbfType;
+   }
+
+   // Convert rbf int to identifier
+   static char GetChar(const int rbfType)
+   {
+      static const char ident[] = { 'G', 'M', 'I',
+                                    'T', 'C',
+                                    '1', '3', '6' };
+      return ident[Check(rbfType)];
+   }
+   
+   // Convert identifier to rbf int
+   static int GetType(const char rbfIdent)
+   {
+      switch (rbfIdent)
+      {
+      case 'G': return Gaussian;
+      case 'M': return Multiquadric;
+      case 'I': return InvMultiquadric;
+      case 'T': return TruncatedGaussian;
+      case 'C': return CompactGaussian;
+      case '1': return Wendland11;
+      case '3': return Wendland31;
+      case '6': return Wendland33;
+      }
+      MFEM_ABORT("unknown RBF identifier: " << rbfIdent);
+      return -1;
+   }
+};
+
+class DistanceMetric
+{
+protected:
+   int Dim;
+public:
+   DistanceMetric(int D) { Dim = D; }
+   virtual ~DistanceMetric() { }
+
+   virtual void SetDim(int D) { Dim = D; }
+   
+   virtual void Distance(const Vector &x,
+                         double &r) const = 0;
+   virtual void DDistance(const Vector &x,
+                          Vector &dr) const = 0;
+   virtual void DDDistance(const Vector &x,
+                           DenseMatrix &ddr) const = 0;
+   
+   static DistanceMetric *GetDistance(int Dim, int pnorm);
+};
+
+class L1Distance : public DistanceMetric
+{
+public:
+   L1Distance(int D) : DistanceMetric(D) { };
+   virtual ~L1Distance() { }
+
+   virtual void Distance(const Vector &x,
+                         double &r) const;
+   virtual void DDistance(const Vector &x,
+                          Vector &dr) const;
+   virtual void DDDistance(const Vector &x,
+                           DenseMatrix &ddr) const;
+};
+
+class L2Distance : public DistanceMetric
+{
+public:
+   L2Distance(int D) : DistanceMetric(D) { };
+   virtual ~L2Distance() { }
+
+   virtual void Distance(const Vector &x,
+                         double &r) const;
+   virtual void DDistance(const Vector &x,
+                          Vector &dr) const;
+   virtual void DDDistance(const Vector &x,
+                           DenseMatrix &ddr) const;
+};
+
+class LpDistance : public DistanceMetric
+{
+   const int p;
+   const double pinv;
+public:
+   LpDistance(int D, int pnorm)
+      : DistanceMetric(D),
+        p(pnorm),
+        pinv(1. / static_cast<double>(p))
+   { };
+   virtual ~LpDistance() { }
+   
+   virtual void Distance(const Vector &x,
+                         double &r) const;
+   virtual void DDistance(const Vector &x,
+                          Vector &dr) const;
+   virtual void DDDistance(const Vector &x,
+                           DenseMatrix &ddr) const;
+};
+
+class KernelFiniteElement : public ScalarFiniteElement
+{
+public:
+   KernelFiniteElement(int D, Geometry::Type G, int Do, int O, int F)
+      : ScalarFiniteElement(D, G, Do, O, F) { }
+   virtual ~KernelFiniteElement() { }
+
+   // Converts integration rule to vector
+   virtual void IntRuleToVec(const IntegrationPoint &ip,
+                             Vector &vec) const;
+
+   virtual bool IsCompact() const = 0;
+   virtual const RBFFunction *Kernel() const = 0;
+
+   // Get range of i,j,k indices that are nonzero for compact support
+   virtual bool TensorIndexed() const { return false; }
+   virtual void GetTensorIndices(const Vector &ip,
+                                 int (&indices)[3][2]) const
+   { MFEM_ABORT("GetTensorIndices(...)"); }
+   virtual void GetTensorNumPoints(int (&tNumPoints)[3]) const
+   { MFEM_ABORT("GetTensorNumPoints(...)"); }
+   
+   using FiniteElement::Project;
+   
+   virtual void Project(Coefficient &coeff,
+                        ElementTransformation &Trans, Vector &dofs) const;
+
+   virtual void Project(const FiniteElement &fe, ElementTransformation &Trans,
+                        DenseMatrix &I) const;
+
+   virtual void GetLocalInterpolation(ElementTransformation &Trans,
+                                      DenseMatrix &I) const
+   { ScalarLocalInterpolation(Trans, I, *this); }
+
+   virtual void GetTransferMatrix(const FiniteElement &fe,
+                                  ElementTransformation &Trans,
+                                  DenseMatrix &I) const
+   { CheckScalarFE(fe).ScalarLocalInterpolation(Trans, I, *this); }
+};
+
+class RBFFiniteElement : public KernelFiniteElement
+{
+private:
+#ifndef MFEM_THREAD_SAFE
+   mutable double r_scr, f_scr, df_scr, ddf_scr;
+   mutable Vector x_scr, y_scr, dy_scr, dr_scr;
+   mutable DenseMatrix ddr_scr;
+   mutable int cInd[3][2];
+#endif
+   bool isCompact; // Is the RBF with the given h compact?
+   int dimPoints[3];
+   int numPointsD; // Number of points across the element in each D
+   double delta; // Distance between points
+   double h; // Shape parameter, approx number of points in 1d support radius
+   double hPhys; // Shape parameter times distance between points times HNorm
+   double hPhysInv; // Inverse hPhys
+   double radPhys; // Radius adjusted by h
+   const RBFFunction *rbf;
+   const DistanceMetric *distance;
+   void InitializeGeometry();
+   
+   virtual void DistanceVec(const int i,
+                            const Vector &x,
+                            Vector &y) const;
+   
+public:
+   RBFFiniteElement(const int D,
+                    const int numPointsD,
+                    const double h,
+                    const int rbfType,
+                    const int distNorm,
+                    const int intOrder);
+   virtual ~RBFFiniteElement() { delete rbf; delete distance; }
+
+   virtual bool TensorIndexed() const { return true; }
+   virtual void GetCompactIndices(const Vector &ip,
+                                  int (&indices)[3][2]) const;
+   virtual void GetGlobalIndices(const Vector &ip,
+                                 int (&indices)[3][2]) const;
+   virtual void GetTensorIndices(const Vector &ip,
+                                 int (&indices)[3][2]) const;
+   virtual void GetTensorNumPoints(int (&tNumPoints)[3]) const
+   {
+      tNumPoints[0] = dimPoints[0];
+      tNumPoints[1] = dimPoints[1];
+      tNumPoints[2] = dimPoints[2];
+   }
+   
+   virtual bool IsCompact() const { return isCompact; }
+   virtual const RBFFunction *Kernel() const { return rbf; }
+   
+   virtual void CalcShape(const IntegrationPoint &ip,
+                          Vector &shape) const;
+   virtual void CalcDShape(const IntegrationPoint &ip,
+                           DenseMatrix &dshape) const;
+   virtual void CalcHessian(const IntegrationPoint &ip,
+                            DenseMatrix &h) const;
+};
+
+class RKFiniteElement : public KernelFiniteElement
+{
+private:
+#ifndef MFEM_THREAD_SAFE
+   mutable double f_scr;
+   mutable Vector x_scr, y_scr, g_scr, c_scr, s_scr, p_scr, df_scr;
+   mutable DenseMatrix q_scr, dq_scr, M_scr;
+   mutable Vector dc_scr[3], dp_scr[3];
+   mutable DenseMatrix dM_scr[3];
+   mutable DenseMatrixInverse Minv_scr;
+   mutable int cInd[3][2];
+   mutable int dimPoints[3];
+#endif
+   int polyOrd, numPoly, numPoly1d;
+   KernelFiniteElement *baseFE;
+
+   virtual void GetPoly(const Vector &x,
+                        Vector &p) const;
+   virtual void GetDPoly(const Vector &x,
+                         Vector &p,
+                         Vector (&dp)[3]) const;
+   virtual void GetG(Vector &g) const;
+   virtual void GetM(const Vector &baseShape,
+                     const IntegrationPoint &ip,
+                     DenseMatrix &M) const;
+   virtual void GetDM(const Vector &baseShape,
+                      const DenseMatrix &baseDeriv,
+                      const IntegrationPoint &ip,
+                      DenseMatrix &M,
+                      DenseMatrix (&dM)[3]) const;
+   virtual void AddToM(const Vector &p,
+                       const double &f,
+                       DenseMatrix &M) const;
+   virtual void AddToDM(const Vector &p,
+                        const Vector (&dp)[3],
+                        const double &f,
+                        const Vector &df,
+                        DenseMatrix (&dM)[3]) const;
+   virtual void CalculateValues(const Vector &c,
+                                const Vector &baseShape,
+                                const IntegrationPoint &ip,
+                                Vector &shape) const;
+   virtual void CalculateDValues(const Vector &c,
+                                 const Vector (&dc)[3],
+                                 const Vector &baseShape,
+                                 const DenseMatrix &baseDShape,
+                                 const IntegrationPoint &ip,
+                                 DenseMatrix &dshape) const;
+   virtual void DistanceVec(const int i,
+                            const Vector &x,
+                            Vector &y) const;
+   
+public:
+   RKFiniteElement(const int D,
+                   const int numPointsD,
+                   const double h,
+                   const int rbfType,
+                   const int distNorm,
+                   const int order,
+                   const int intOrder);
+   virtual ~RKFiniteElement() { delete baseFE; }
+   
+   virtual bool IsCompact() const { return baseFE->IsCompact(); }
+   virtual const RBFFunction *Kernel() const { return baseFE->Kernel(); }
+   
+   static int GetNumPoly(int polyOrd, int dim);
+   
+   virtual void CalcShape(const IntegrationPoint &ip,
+                          Vector &shape) const;
+   virtual void CalcDShape(const IntegrationPoint &ip,
+                           DenseMatrix &dshape) const;
+   // Should put in a method to calculate shape and dshape simultaneously
+};
+
 } // namespace mfem
 
 #endif

fem/fe_coll.cpp

@@ -278,6 +278,35 @@ FiniteElementCollection *FiniteElementCollection::New(const char *name)
          fec = new NURBSFECollection();
       }
    }
+   else if (!strncmp(name, "RBF", 3) || !strncmp(name, "RK", 2))
+   {
+      // Example: RK4_G_2_V_2D_0020_4.01
+      // (RK order 4, Gaussian, L2 dist, Value map, 2 dimensions,
+      //  20 points across element, smoothing length of 4.01)
+      const int dim = atoi(name + 10);
+      const int numPoints = atoi(name + 13);
+      const double h = atof(name + 18);
+      const int rbfType = RBFType::GetType(name[4]);
+      const int distNorm = atoi(name + 6);
+      const int mapType = (name[8] == 'V'
+                           ? FiniteElement::VALUE
+                           : FiniteElement::INTEGRAL);
+      const int intOrder = 2; // fix this for now
+      
+      if (!strncmp(name, "RK", 2))
+      {
+         int order = atoi(name + 2);
+         fec = new KernelFECollection(dim, numPoints, h,
+                                      rbfType, distNorm, order,
+                                      intOrder, mapType);
+      }
+      else
+      {
+         fec = new KernelFECollection(dim, numPoints, h,
+                                      rbfType, distNorm, -1,
+                                      intOrder, mapType);
+      }
+   }
    else
    {
       MFEM_ABORT("unknown FiniteElementCollection: " << name);
@@ -2618,4 +2647,175 @@ FiniteElementCollection *NURBSFECollection::GetTraceCollection() const
    return NULL;
 }
 
+KernelFECollection::KernelFECollection(const int dim,
+                                       const int numPointsD,
+                                       const double h,
+                                       const int rbfType,
+                                       const int distNorm,
+                                       const int order,
+                                       const int intOrder,
+                                       const int mapType)
+{
+   const char *mapStr = NULL;
+   switch (mapType)
+   {
+      case FiniteElement::VALUE:    mapStr = "V"; break;
+      case FiniteElement::INTEGRAL: mapStr = "I"; break;
+      default:
+         MFEM_ABORT("invalid mapType: " << mapType);
+   }
+   if (order == -1)
+   {
+      snprintf(d_name, 32, "RBF_%c_%d_%s_%dD_%04d_%.2f",
+               (int)RBFType::GetChar(rbfType), distNorm,
+               mapStr, dim, numPointsD, h);
+   }
+   else if (order >= 0)
+   {
+      snprintf(d_name, 32, "RK%d_%c_%d_%s_%dD_%04d_%.2f", order,
+               (int)RBFType::GetChar(rbfType), distNorm,
+               mapStr, dim, numPointsD, h);
+   }
+   else
+   {
+      MFEM_ABORT("invalid order: " << order);
+   }
+   
+   for (int g = 0; g < Geometry::NumGeom; ++g)
+   {
+      L2_Elements[g] = NULL;
+      Tr_Elements[g] = NULL;
+   }
+   for (int i = 0; i < 2; i++)
+   {
+      SegDofOrd[i] = NULL;
+   }
+   OtherDofOrd = NULL;
+   
+   if (dim == 0)
+   {
+      L2_Elements[Geometry::POINT] = new PointFiniteElement;
+   }
+   else if (dim == 1)
+   {
+      if (order == -1)
+      {
+         L2_Elements[Geometry::SEGMENT]
+            = new RBFFiniteElement(1, numPointsD, h,
+                                   rbfType, distNorm, intOrder);
+      }
+      else {
+         L2_Elements[Geometry::SEGMENT]
+            = new RKFiniteElement(1, numPointsD, h,
+                                  rbfType, distNorm, order, intOrder);
+      }
+      L2_Elements[Geometry::SEGMENT]->SetMapType(mapType);
+      Tr_Elements[Geometry::POINT] = new PointFiniteElement;
+   }
+   else if (dim == 2)
+   {
+      if (order == -1)
+      {
+         L2_Elements[Geometry::SQUARE]
+            = new RBFFiniteElement(2, numPointsD, h,
+                                   rbfType, distNorm, intOrder);
+         Tr_Elements[Geometry::SEGMENT]
+            = new RBFFiniteElement(1, numPointsD, h,
+                                   rbfType, distNorm, intOrder);
+      }
+      else {
+         L2_Elements[Geometry::SQUARE]
+            = new RKFiniteElement(2, numPointsD, h,
+                                  rbfType, distNorm, order, intOrder);
+         Tr_Elements[Geometry::SEGMENT]
+            = new RKFiniteElement(1, numPointsD, h,
+                                  rbfType, distNorm, order, intOrder);
+      }
+      L2_Elements[Geometry::SQUARE]->SetMapType(mapType);
+   }
+   else if (dim == 3)
+   {
+      if (order == -1)
+      {
+         L2_Elements[Geometry::CUBE]
+            = new RBFFiniteElement(3, numPointsD, h,
+                                   rbfType, distNorm, intOrder);
+         Tr_Elements[Geometry::SQUARE]
+            = new RBFFiniteElement(2, numPointsD, h,
+                                   rbfType, distNorm, intOrder);
+      }
+      else {
+         L2_Elements[Geometry::CUBE]
+            = new RKFiniteElement(3, numPointsD, h,
+                                  rbfType, distNorm, order, intOrder);
+         Tr_Elements[Geometry::SQUARE]
+            = new RKFiniteElement(2, numPointsD, h,
+                                  rbfType, distNorm, order, intOrder);
+      }
+      L2_Elements[Geometry::CUBE]->SetMapType(mapType);
+   }
+   
+   if (dim == 1)
+   {
+      SegDofOrd[0] = new int[2*numPointsD];
+      SegDofOrd[1] = SegDofOrd[0] + numPointsD;
+      for (int i = 0; i < numPointsD; ++i) {
+         SegDofOrd[0][i] = i;
+         SegDofOrd[1][i] = numPointsD - i - 1;
+      }
+   }
+   else
+   {
+      const int geomType = TensorBasisElement::GetTensorProductGeometry(dim);
+      const int dof = L2_Elements[geomType]->GetDof();
+      OtherDofOrd = new int[dof];
+      for (int i = 0; i < dof; ++i)
+      {
+         OtherDofOrd[i] = i;
+      }
+   }
 }
+
+KernelFECollection::~KernelFECollection()
+{
+   delete [] OtherDofOrd;
+   delete [] SegDofOrd[0];
+   for (int i = 0; i < Geometry::NumGeom; ++i)
+   {
+      delete L2_Elements[i];
+   }
+}
+
+const FiniteElement *
+KernelFECollection::FiniteElementForGeometry(Geometry::Type GeomType) const
+{
+   return L2_Elements[GeomType];
+}
+
+const FiniteElement *
+KernelFECollection::TraceFiniteElementForGeometry(Geometry::Type GeomType) const
+{
+   return Tr_Elements[GeomType];
+}
+
+int KernelFECollection::DofForGeometry(Geometry::Type GeomType) const
+{
+   if (L2_Elements[GeomType])
+   {
+      return L2_Elements[GeomType]->GetDof();
+   }
+   return 0;
+}
+
+const int *KernelFECollection::DofOrderForOrientation(Geometry::Type GeomType, int Or) const
+{
+   if (GeomType == Geometry::SEGMENT) {
+      return (Or > 0) ? SegDofOrd[0] : SegDofOrd[1];
+   }
+   else
+   {
+      return (Or == 0) ? OtherDofOrd : NULL;
+   }
+}
+
+} // namespace mfem

fem/fe_coll.hpp

@@ -895,6 +895,43 @@ public:
    virtual ~Local_FECollection() { delete Local_Element; }
 };
 
-}
+/// Radial basis function collection
+class KernelFECollection : public FiniteElementCollection
+{
+private:
+   int maxDim;
+   char d_name[32];
+   ScalarFiniteElement *Tr_Elements[Geometry::NumGeom];
+   ScalarFiniteElement *L2_Elements[Geometry::NumGeom];
+   int *SegDofOrd[2]; // for rotating segment dofs in 1D
+   int *OtherDofOrd; 
+   bool ValidGeomType(Geometry::Type GeomType) const;
+   
+public:
+   KernelFECollection(const int D,
+                      const int numPointsD,
+                      const double h,
+                      const int rbfType,
+                      const int distNorm,
+                      const int order = -1,
+                      const int intOrder = 2, // num integration points per 1d point
+                      const int mapType = FiniteElement::VALUE);
+   virtual ~KernelFECollection();
+   
+   virtual const FiniteElement *
+   FiniteElementForGeometry(Geometry::Type GeomType) const;
+
+   virtual const FiniteElement *
+   TraceFiniteElementForGeometry(Geometry::Type GeomType) const;
+
+   virtual int DofForGeometry(Geometry::Type GeomType) const;
+
+   virtual const int * DofOrderForOrientation(Geometry::Type GeomType, int Or) const;
+
+   virtual const char * Name() const { return d_name; }
+   // virtual int GetContType() const { return DISCONTINUOUS; }
+};
+
+} // namespace mfem
 
 #endif

fem/lininteg.cpp

@@ -64,6 +64,57 @@ void DomainLFIntegrator::AssembleDeltaElementVect(
 }
 
 
+void DomainLFGradIntegrator::AssembleRHSElementVect(
+    const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
+{
+   int dim = el.GetDim();
+   int dof = el.GetDof();
+   
+   Vector Qvec(dim);
+   dshape.SetSize(dof, dim);
+   dshapeQ.SetSize(dof);
+   elvect.SetSize(dof);
+   elvect = 0.0;
+
+
+   const IntegrationRule *ir = IntRule;
+   if (ir == NULL)
+   {
+      ir = &IntRules.Get(el.GetGeomType(), oa * el.GetOrder() + ob);
+   }
+   
+   for (int i = 0; i < ir->GetNPoints(); i++)
+   {
+      const IntegrationPoint &ip = ir->IntPoint(i);
+
+      Tr.SetIntPoint(&ip);
+      Q.Eval(Qvec, Tr, ip);
+      el.CalcDShape(ip, dshape);
+      
+      dshape.Mult(Qvec, dshapeQ);
+      elvect.Add(ip.weight * Tr.Weight(), dshapeQ);
+   }
+}
+
+void DomainLFGradIntegrator::AssembleDeltaElementVect(
+   const FiniteElement &fe, ElementTransformation &Trans, Vector &elvect)
+{
+   MFEM_ASSERT(delta != NULL, "coefficient must be DeltaCoefficient");
+   
+   int dim = fe.GetDim();
+   int dof = fe.GetDof();
+   
+   Vector Qvec(dim);
+   elvect.SetSize(dof);
+   dshape.SetSize(dof, dim);
+   
+   fe.CalcPhysDShape(Trans, dshape);
+   vec_delta->EvalDelta(Qvec, Trans, Trans.GetIntPoint());
+   
+   dshape.Mult(Qvec, elvect);
+}
+
+
 void BoundaryLFIntegrator::AssembleRHSElementVect(
    const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
 {

fem/lininteg.hpp

@@ -119,6 +119,39 @@ public:
    using LinearFormIntegrator::AssembleRHSElementVect;
 };
 
+/// Class for domain integration L(v) := (f, grad v)
+class DomainLFGradIntegrator : public DeltaLFIntegrator
+{
+   DenseMatrix dshape;
+   Vector dshapeQ;
+   VectorCoefficient &Q;
+   int oa, ob;
+public:
+   /// Constructs a domain integrator with a given Coefficient
+   DomainLFGradIntegrator(VectorCoefficient &QF, int a = 2, int b = 0)
+   // the old default was a = 1, b = 1
+   // for simple elliptic problems a = 2, b = -2 is OK
+      : DeltaLFIntegrator(QF), Q(QF), oa(a), ob(b) { }
+
+   /// Constructs a domain integrator with a given Coefficient
+   DomainLFGradIntegrator(VectorCoefficient &QF, const IntegrationRule *ir)
+      : DeltaLFIntegrator(QF, ir), Q(QF), oa(1), ob(1) { }
+
+   /** Given a particular Finite Element and a transformation (Tr)
+       computes the element right hand side element vector, elvect. */
+   virtual void AssembleRHSElementVect(const FiniteElement &el,
+                                       ElementTransformation &Tr,
+                                       Vector &elvect);
+
+   virtual void AssembleDeltaElementVect(const FiniteElement &fe,
+                                         ElementTransformation &Trans,
+                                         Vector &elvect);
+
+   void ResetCoefficient(VectorCoefficient &q) { Q = q;  }
+
+   using LinearFormIntegrator::AssembleRHSElementVect;
+};
+
 /// Class for boundary integration L(v) := (g, v)
 class BoundaryLFIntegrator : public LinearFormIntegrator
 {

linalg/blockoperator.cpp

@@ -326,4 +326,207 @@ BlockLowerTriangularPreconditioner::~BlockLowerTriangularPreconditioner()
    }
 }
 
+BlockUpperTriangularPreconditioner::BlockUpperTriangularPreconditioner(
+   const Array<int> & offsets_)
+   : Solver(offsets_.Last()),
+     owns_blocks(0),
+     nBlocks(offsets_.Size() - 1),
+     offsets(0),
+     op(nBlocks, nBlocks)
+{
+   op = static_cast<Operator *>(NULL);
+   offsets.MakeRef(offsets_);
+}
+
+void BlockUpperTriangularPreconditioner::SetDiagonalBlock(int iblock,
+                                                          Operator *op)
+{
+   MFEM_VERIFY(offsets[iblock+1] - offsets[iblock] == op->Height() &&
+               offsets[iblock+1] - offsets[iblock] == op->Width(),
+               "incompatible Operator dimensions");
+
+   SetBlock(iblock, iblock, op);
+}
+
+void BlockUpperTriangularPreconditioner::SetBlock(int iRow, int iCol,
+                                                  Operator *opt)
+{
+   MFEM_VERIFY(iRow <= iCol,"cannot set block in lower triangle");
+   MFEM_VERIFY(offsets[iRow+1] - offsets[iRow] == opt->NumRows() &&
+               offsets[iCol+1] - offsets[iCol] == opt->NumCols(),
+               "incompatible Operator dimensions");
+
+   op(iRow, iCol) = opt;
+}
+
+// Operator application
+void BlockUpperTriangularPreconditioner::MultTranspose (const Vector & x,
+                                               Vector & y) const
+{
+   MFEM_ASSERT(x.Size() == width, "incorrect input Vector size");
+   MFEM_ASSERT(y.Size() == height, "incorrect output Vector size");
+
+   yblock.Update(y.GetData(),offsets);
+   xblock.Update(x.GetData(),offsets);
+
+   y = 0.0;
+   for (int iRow=0; iRow < nBlocks; ++iRow)
+   {
+      tmp.SetSize(offsets[iRow+1] - offsets[iRow]);
+      tmp2.SetSize(offsets[iRow+1] - offsets[iRow]);
+      tmp2 = 0.0;
+      tmp2 += xblock.GetBlock(iRow);
+      for (int jCol=0; jCol < iRow; ++jCol)
+      {
+         if (op(iRow,jCol))
+         {
+            op(iRow,jCol)->MultTranspose(yblock.GetBlock(jCol), tmp);
+            tmp2 -= tmp;
+         }
+      }
+      if (op(iRow,iRow))
+      {
+         op(iRow,iRow)->MultTranspose(tmp2, yblock.GetBlock(iRow));
+      }
+      else
+      {
+         yblock.GetBlock(iRow) = tmp2;
+      }
+   }
+}
+
+// Action of the transpose operator
+void BlockUpperTriangularPreconditioner::Mult(const Vector & x,
+                                                        Vector & y) const
+{
+   MFEM_ASSERT(x.Size() == height, "incorrect input Vector size");
+   MFEM_ASSERT(y.Size() == width, "incorrect output Vector size");
+
+   yblock.Update(y.GetData(),offsets);
+   xblock.Update(x.GetData(),offsets);
+
+   y = 0.0;
+   for (int iRow=nBlocks-1; iRow >=0; --iRow)
+   {
+      tmp.SetSize(offsets[iRow+1] - offsets[iRow]);
+      tmp2.SetSize(offsets[iRow+1] - offsets[iRow]);
+      tmp2 = 0.0;
+      tmp2 += xblock.GetBlock(iRow);
+      for (int jCol=iRow+1; jCol < nBlocks; ++jCol)
+      {
+         if (op(jCol,iRow))
+         {
+            op(jCol,iRow)->Mult(yblock.GetBlock(jCol), tmp);
+            tmp2 -= tmp;
+         }
+      }
+      if (op(iRow,iRow))
+      {
+         op(iRow,iRow)->Mult(tmp2, yblock.GetBlock(iRow));
+      }
+      else
+      {
+         yblock.GetBlock(iRow) = tmp2;
+      }
+   }
+}
+
+BlockUpperTriangularPreconditioner::~BlockUpperTriangularPreconditioner()
+{
+   if (owns_blocks)
+   {
+      for (int iRow=0; iRow < nBlocks; ++iRow)
+      {
+         for (int jCol=0; jCol < nBlocks; ++jCol)
+         {
+            delete op(jCol,iRow);
+         }
+      }
+   }
+}
+
+BlockLDUPreconditioner::BlockLDUPreconditioner(
+   const Array<int> & offsets_, int schur_index_)
+   : Solver(offsets_.Last()),
+     owns_blocks(0),
+     nBlocks(offsets_.Size() - 1),
+     offsets(0),
+     op(nBlocks, nBlocks),
+     schur_index(schur_index_)
+{
+   op = static_cast<Operator *>(NULL);
+   offsets.MakeRef(offsets_);
+}
+
+void BlockLDUPreconditioner::SetDiagonalBlock(int iblock, Operator *op)
+{
+   MFEM_VERIFY(offsets[iblock+1] - offsets[iblock] == op->Height() &&
+               offsets[iblock+1] - offsets[iblock] == op->Width(),
+               "incompatible Operator dimensions");
+   MFEM_VERIFY(iblock >=0 && iblock <= 1,
+               "Only valid for 2x2 block matrices");
+
+   SetBlock(iblock, iblock, op);
+}
+
+void BlockLDUPreconditioner::SetBlock(int iRow, int iCol, Operator *opt)
+{
+   MFEM_VERIFY(iRow >=0 && iRow <= 1, "Only valid for 2x2 block matrices");
+   MFEM_VERIFY(iCol >=0 && iCol <= 1, "Only valid for 2x2 block matrices");
+   MFEM_VERIFY(offsets[iRow+1] - offsets[iRow] == opt->NumRows() &&
+               offsets[iCol+1] - offsets[iCol] == opt->NumCols(),
+               "incompatible Operator dimensions");
+
+   op(iRow, iCol) = opt;
+}
+
+// Action of the transpose operator
+void BlockLDUPreconditioner::Mult(const Vector & x, Vector & y) const
+{
+   MFEM_ASSERT(x.Size() == height, "incorrect input Vector size");
+   MFEM_ASSERT(y.Size() == width, "incorrect output Vector size");
+
+   yblock.Update(y.GetData(),offsets);
+   xblock.Update(x.GetData(),offsets);
+   y = 0.0;
+
+   // Schur-complement in the (1,1) block (i.e., (0,0) block w/ zero indexing)
+   // TODO
+
+   // Schur-complement in the (2,2) block (i.e., (1,1) block w/ zero indexing)
+   // A^{-1} = [I, -A11^{-1}A12; 0, I] * [I, 0; 0, S22^{-2}] * 
+   //          [I, 0; -A21, I] * [A11^{-1}, 0; 0, I]
+   op(0,0) -> Mult(xblock.GetBlock(0), yblock.GetBlock(0));
+   
+   tmp.SetSize(offsets[2] - offsets[1]);
+   op(1,0) -> Mult(yblock.GetBlock(0), tmp);
+   tmp *= -1;
+   tmp += xblock.GetBlock(1);
+
+   op(1,1) -> Mult(tmp, yblock.GetBlock(1));
+
+   // BUG IN THIS SECTION
+   tmp3.SetSize(offsets[1] - offsets[0]);
+   tmp2.SetSize(offsets[1] - offsets[0]);
+   tmp2 = 0.0;
+   op(0,1) -> Mult(yblock.GetBlock(1), tmp3);
+   op(0,0) -> Mult(tmp3, tmp2);
+   tmp = yblock.GetBlock(0);
+   yblock.GetBlock(0) = tmp - tmp2; // <--- this line
+}
+
+BlockLDUPreconditioner::~BlockLDUPreconditioner()
+{
+   if (owns_blocks)
+   {
+      for (int iRow=0; iRow < nBlocks; ++iRow)
+      {
+         for (int jCol=0; jCol < nBlocks; ++jCol)
+         {
+            delete op(jCol,iRow);
+         }
+      }
+   }
+}
+
 }

linalg/blockoperator.hpp

@@ -212,7 +212,7 @@ public:
     *  @note BlockLowerTriangularPreconditioner will not own/copy the data
     *  contained in @a offsets.
     */
-   BlockLowerTriangularPreconditioner(const Array<int> & offsets);
+   BlockLowerTriangularPreconditioner(const Array<int> & offsets_);
 
    //! Add block op in the block-entry (iblock, iblock).
    /**
@@ -268,5 +268,179 @@ private:
    mutable Vector tmp2;
 };
 
+//! @class BlockUpperTriangularPreconditioner
+/**
+ * \brief A class to handle Block upper triangular preconditioners in a
+ * matrix-free implementation.
+ *
+ * Usage:
+ * - Use the constructors to define the block structure
+ * - Use SetBlock() to fill the BlockOperator
+ * - Diagonal blocks of the preconditioner should approximate the inverses of
+ *   the diagonal block of the matrix
+ * - Off-diagonal blocks of the preconditioner should match/approximate those of
+ *   the original matrix
+ * - Use the method Mult() and MultTranspose() to apply the operator to a vector.
+ *
+ * If a diagonal block is not set, it is assumed to be an identity block, if an
+ * off-diagonal block is not set, it is assumed to be a zero block.
+ *
+ */
+class BlockUpperTriangularPreconditioner : public Solver
+{
+public:
+   //! Constructor for BlockUpperTriangularPreconditioners with the same
+   //! block-structure for rows and columns.
+   /**
+    *  @param offsets  Offsets that mark the start of each row/column block
+    *                  (size nBlocks+1).
+    *
+    *  @note BlockUpperTriangularPreconditioner will not own/copy the data
+    *  contained in @a offsets.
+    */
+   BlockUpperTriangularPreconditioner(const Array<int> & offsets_);
+
+   //! Add block op in the block-entry (iblock, iblock).
+   /**
+    * @param iblock  The block will be inserted in location (iblock, iblock).
+    * @param op      The Operator to be inserted.
+    */
+   void SetDiagonalBlock(int iblock, Operator *op);
+   //! Add a block op in the block-entry (iblock, jblock).
+   /**
+    * @param iRow, iCol  The block will be inserted in location (iRow, iCol).
+    * @param op          The Operator to be inserted.
+    */
+   void SetBlock(int iRow, int iCol, Operator *op);
+   //! This method is present since required by the abstract base class Solver
+   virtual void SetOperator(const Operator &op) { }
+
+   //! Return the number of blocks
+   int NumBlocks() const { return nBlocks; }
+
+   //! Return a reference to block i,j.
+   Operator & GetBlock(int iblock, int jblock)
+   { MFEM_VERIFY(op(iblock,jblock), ""); return *op(iblock,jblock); }
+
+   //! Return the offsets for block starts
+   Array<int> & Offsets() { return offsets; }
+
+   /// Operator application
+   virtual void Mult (const Vector & x, Vector & y) const;
+
+   /// Action of the transpose operator
+   virtual void MultTranspose (const Vector & x, Vector & y) const;
+
+   ~BlockUpperTriangularPreconditioner();
+
+   //! Controls the ownership of the blocks: if nonzero,
+   //! BlockUpperTriangularPreconditioner will delete all blocks that are set
+   //! (non-NULL); the default value is zero.
+   int owns_blocks;
+
+private:
+   //! Number of block rows/columns
+   int nBlocks;
+   //! Offsets for the starting position of each block
+   Array<int> offsets;
+   //! 2D array that stores each block of the operator.
+   Array2D<Operator *> op;
+
+   //! Temporary Vectors used to efficiently apply the Mult and MultTranspose
+   //! methods.
+   mutable BlockVector xblock;
+   mutable BlockVector yblock;
+   mutable Vector tmp;
+   mutable Vector tmp2;
+};
+
+//! @class BlockLDUPreconditioner
+/**
+ * \brief A class to handle Block upper triangular preconditioners in a
+ * matrix-free implementation.
+ *
+ * Usage:
+ * - Use the constructors to define the block structure
+ * - Use SetBlock() to fill the BlockOperator
+ * - Diagonal blocks of the preconditioner should approximate the inverses of
+ *   the diagonal block of the matrix
+ * - Off-diagonal blocks of the preconditioner should match/approximate those of
+ *   the original matrix
+ * - Use the method Mult() and MultTranspose() to apply the operator to a vector.
+ *
+ * If a diagonal block is not set, it is assumed to be an identity block, if an
+ * off-diagonal block is not set, it is assumed to be a zero block.
+ *
+ */
+class BlockLDUPreconditioner : public Solver
+{
+public:
+   //! Constructor for BlockLDUPreconditioners with the same
+   //! block-structure for rows and columns.
+   /**
+    *  @param offsets  Offsets that mark the start of each row/column block
+    *                  (size nBlocks+1).
+    *
+    *  @note BlockLDUPreconditioner will not own/copy the data
+    *  contained in @a offsets.
+    */
+   BlockLDUPreconditioner(const Array<int> & offsets_, int schur_index_=1);
+
+   //! Add block op in the block-entry (iblock, iblock).
+   /**
+    * @param iblock  The block will be inserted in location (iblock, iblock).
+    * @param op      The Operator to be inserted.
+    */
+   void SetDiagonalBlock(int iblock, Operator *op);
+   //! Add a block op in the block-entry (iblock, jblock).
+   /**
+    * @param iRow, iCol  The block will be inserted in location (iRow, iCol).
+    * @param op          The Operator to be inserted.
+    */
+   void SetBlock(int iRow, int iCol, Operator *op);
+   //! This method is present since required by the abstract base class Solver
+   virtual void SetOperator(const Operator &op) { }
+
+   //! Return the number of blocks
+   int NumBlocks() const { return nBlocks; }
+
+   //! Return a reference to block i,j.
+   Operator & GetBlock(int iblock, int jblock)
+   { MFEM_VERIFY(op(iblock,jblock), ""); return *op(iblock,jblock); }
+
+   //! Return the offsets for block starts
+   Array<int> & Offsets() { return offsets; }
+
+   /// Operator application
+   virtual void Mult (const Vector & x, Vector & y) const;
+
+   /// Action of the transpose operator
+   virtual void MultTranspose (const Vector & x, Vector & y) const
+   { MFEM_WARNING("MultTransport not implemented for LDU.\n;"); }
+
+   ~BlockLDUPreconditioner();
+
+   //! Controls the ownership of the blocks: if nonzero,
+   //! BlockLDUPreconditioner will delete all blocks that are set
+   //! (non-NULL); the default value is zero.
+   int owns_blocks;
+
+private:
+   //! Number of block rows/columns
+   int nBlocks, schur_index;
+   //! Offsets for the starting position of each block
+   Array<int> offsets;
+   //! 2D array that stores each block of the operator.
+   Array2D<Operator *> op;
+
+   //! Temporary Vectors used to efficiently apply the Mult and MultTranspose
+   //! methods.
+   mutable BlockVector xblock;
+   mutable BlockVector yblock;
+   mutable Vector tmp;
+   mutable Vector tmp2;
+   mutable Vector tmp3;
+};
+
 }
 #endif /* MFEM_BLOCKOPERATOR */

linalg/hypre.cpp

@@ -173,6 +173,13 @@ HypreParVector::HypreParVector(ParFiniteElementSpace *pfes)
    own_ParVector = 1;
 }
 
+void HypreParVector::WrapHypreParVector(hypre_ParVector *y)
+{
+   x = y;
+   _SetDataAndSize_();
+   own_ParVector = 0;
+}
+
 Vector * HypreParVector::GlobalVector() const
 {
    hypre_Vector *hv = hypre_ParVectorToVectorAll(*this);
@@ -967,17 +974,26 @@ static void MakeWrapper(const hypre_CSRMatrix *mat, SparseMatrix &wrapper)
    wrapper.Swap(tmp);
 }
 
+
 void HypreParMatrix::GetDiag(SparseMatrix &diag) const
 {
    MakeWrapper(A->diag, diag);
 }
 
+
 void HypreParMatrix::GetOffd(SparseMatrix &offd, HYPRE_Int* &cmap) const
 {
    MakeWrapper(A->offd, offd);
    cmap = A->col_map_offd;
 }
 
+
+void HypreParMatrix::GetProcRows(SparseMatrix &colCSRMat)
+{
+   MakeWrapper(hypre_MergeDiagAndOffd(A), colCSRMat);
+}
+
+
 void HypreParMatrix::GetBlocks(Array2D<HypreParMatrix*> &blocks,
                                bool interleaved_rows,
                                bool interleaved_cols) const
@@ -1018,6 +1034,46 @@ HypreParMatrix * HypreParMatrix::Transpose() const
    return new HypreParMatrix(At);
 }
 
+HypreParMatrix * HypreParMatrix::ExtractSubmatrix(Array<int> &indices,
+                                                  double threshhold) const
+{
+   if (!(A->comm))
+   {
+      BuildComm();
+   }
+
+   hypre_ParCSRMatrix *submat;
+
+   // Get number of rows stored on this processor
+   int local_num_vars = hypre_CSRMatrixNumRows(hypre_ParCSRMatrixDiag(A));
+
+   // Form hypre CF-splitting array designating submatrix as F-points (-1)
+   int *CF_marker = new int[local_num_vars];
+   std::fill_n(CF_marker, local_num_vars, 1);
+   for (int j=0; j<indices.Size(); j++)
+   {
+      if (indices[j] > local_num_vars)
+      {
+         MFEM_WARNING("WARNING : " << indices[j] << " > " << local_num_vars);
+      }
+      CF_marker[indices[j]] = -1;
+   }
+
+   // Construct cpts_global array on hypre matrix structure
+   int *cpts_global;
+   hypre_BoomerAMGCoarseParms(MPI_COMM_WORLD, local_num_vars, 1, NULL,
+                              CF_marker, NULL, &cpts_global);
+
+   // Extract submatrix into *submat
+   hypre_ParCSRMatrixExtractSubmatrixFC(A, CF_marker, cpts_global,
+                                        "FF", &submat, threshhold);
+
+   delete[] CF_marker;
+   free(cpts_global);
+   return new HypreParMatrix(submat);
+}
+
+
 HYPRE_Int HypreParMatrix::Mult(HypreParVector &x, HypreParVector &y,
                                double a, double b)
 {
@@ -1590,21 +1646,68 @@ void HypreParMatrix::Destroy()
    }
 }
 
+/* job =  0, extract block diagonal of A and scale A into C
+ * job =  1, job 0 + scale b into d
+ * job =  2, use A to scale b only
+ */
+int BlockInvScal(const HypreParMatrix *A, HypreParMatrix *C,
+                 const Vector *b, HypreParVector *d, int block, int job)
+{
+   if (0 == job || 1 == job)
+   {
+      hypre_ParCSRMatrix *C_hypre;
+      hypre_ParcsrBdiagInvScal(*A, block, &C_hypre);
+      /* XXX: FIXME drop in BdiagInvScal */
+      hypre_ParCSRMatrixDropSmallEntries(C_hypre, 1e-15, 1);
+      (*C).WrapHypreParCSRMatrix(C_hypre);
+   }
+
+   if (1 == job || 2 == job)
+   {
+      HypreParVector *b_Hypre = new HypreParVector(A->GetComm(),
+                                                   A->GetGlobalNumRows(),
+                                                   b->GetData(), A->GetRowStarts());
+      hypre_ParVector *d_hypre;
+      hypre_ParvecBdiagInvScal(*b_Hypre, block, &d_hypre, *A);
+      delete b_Hypre;
+
+      d->WrapHypreParVector(d_hypre);
+      d->SetOwnership(true);
+
+      return 0;
+   }
+
+   return -1;
+}
+
 HypreParMatrix *Add(double alpha, const HypreParMatrix &A,
                     double beta,  const HypreParMatrix &B)
 {
-   hypre_ParCSRMatrix *C_hypre =
-      internal::hypre_ParCSRMatrixAdd(const_cast<HypreParMatrix &>(A),
-                                      const_cast<HypreParMatrix &>(B));
-   MFEM_VERIFY(C_hypre, "error in hypre_ParCSRMatrixAdd");
+   hypre_ParCSRMatrix *C;
+   hypre_ParcsrAdd(alpha, A, beta, B, &C);
+   hypre_MatvecCommPkgCreate(C);
 
+   return new HypreParMatrix(C);
+}
+
+HypreParMatrix *HypreParMatrixAdd(double alpha, const HypreParMatrix &A,
+                                  double beta,  const HypreParMatrix &B)
+{
+   hypre_ParCSRMatrix *C_hypre;
+   hypre_ParcsrAdd(alpha, A, beta, B, &C_hypre);
    hypre_MatvecCommPkgCreate(C_hypre);
-   HypreParMatrix *C = new HypreParMatrix(C_hypre);
-   *C = 0.0;
-   C->Add(alpha, A);
-   C->Add(beta, B);
 
-   return C;
+   return new HypreParMatrix(C_hypre);
+}
+
+HypreParMatrix * ParAdd(const HypreParMatrix *A, const HypreParMatrix *B)
+{
+   hypre_ParCSRMatrix *C;
+   hypre_ParcsrAdd(1.0, *A, 1.0, *B, &C);
+
+   hypre_MatvecCommPkgCreate(C);
+
+   return new HypreParMatrix(C);
 }
 
 HypreParMatrix * ParMult(const HypreParMatrix *A, const HypreParMatrix *B,
@@ -1624,15 +1727,6 @@ HypreParMatrix * ParMult(const HypreParMatrix *A, const HypreParMatrix *B,
    return C;
 }
 
-HypreParMatrix * ParAdd(const HypreParMatrix *A, const HypreParMatrix *B)
-{
-   hypre_ParCSRMatrix * C = internal::hypre_ParCSRMatrixAdd(*A,*B);
-
-   hypre_MatvecCommPkgCreate(C);
-
-   return new HypreParMatrix(C);
-}
-
 HypreParMatrix * RAP(const HypreParMatrix *A, const HypreParMatrix *P)
 {
    HYPRE_Int P_owns_its_col_starts =
@@ -1891,9 +1985,9 @@ HypreSmoother::HypreSmoother(HypreParMatrix &_A, int _type,
    SetOperator(_A);
 }
 
-void HypreSmoother::SetType(HypreSmoother::Type _type, int _relax_times)
+void HypreSmoother::SetType(int _type, int _relax_times)
 {
-   type = static_cast<int>(_type);
+   type = _type;
    relax_times = _relax_times;
 }
 
@@ -2176,6 +2270,8 @@ HypreSolver::HypreSolver()
 {
    A = NULL;
    setup_called = 0;
+   final_res_norm = -1;
+   num_iterations = -1;
    B = X = NULL;
    error_mode = ABORT_HYPRE_ERRORS;
 }
@@ -2185,6 +2281,8 @@ HypreSolver::HypreSolver(HypreParMatrix *_A)
 {
    A = _A;
    setup_called = 0;
+   final_res_norm = -1;
+   num_iterations = -1;
    B = X = NULL;
    error_mode = ABORT_HYPRE_ERRORS;
 }
@@ -2349,8 +2447,6 @@ void HyprePCG::Mult(const HypreParVector &b, HypreParVector &x) const
 {
    int myid;
    HYPRE_Int time_index = 0;
-   HYPRE_Int num_iterations;
-   double final_res_norm;
    MPI_Comm comm;
    HYPRE_Int print_level;
 
@@ -2394,6 +2490,9 @@ void HyprePCG::Mult(const HypreParVector &b, HypreParVector &x) const
    x.HostReadWrite();
 
    HYPRE_ParCSRPCGSolve(pcg_solver, *A, b, x);
+   HYPRE_ParCSRPCGGetNumIterations(pcg_solver, &num_iterations);
+   HYPRE_ParCSRPCGGetFinalRelativeResidualNorm(pcg_solver,
+                                               &final_res_norm);
 
    if (print_level > 0)
    {
@@ -2405,10 +2504,6 @@ void HyprePCG::Mult(const HypreParVector &b, HypreParVector &x) const
          hypre_ClearTiming();
       }
 
-      HYPRE_ParCSRPCGGetNumIterations(pcg_solver, &num_iterations);
-      HYPRE_ParCSRPCGGetFinalRelativeResidualNorm(pcg_solver,
-                                                  &final_res_norm);
-
       MPI_Comm_rank(comm, &myid);
 
       if (myid == 0)
@@ -2483,6 +2578,12 @@ void HypreGMRES::SetTol(double tol)
    HYPRE_GMRESSetTol(gmres_solver, tol);
 }
 
+void HypreGMRES::SetAbsTol(double tol)
+{
+   HYPRE_GMRESSetTol(gmres_solver, 0.0);
+   HYPRE_GMRESSetAbsoluteTol(gmres_solver, tol);
+}
+
 void HypreGMRES::SetMaxIter(int max_iter)
 {
    HYPRE_GMRESSetMaxIter(gmres_solver, max_iter);
@@ -2517,8 +2618,6 @@ void HypreGMRES::Mult(const HypreParVector &b, HypreParVector &x) const
 {
    int myid;
    HYPRE_Int time_index = 0;
-   HYPRE_Int num_iterations;
-   double final_res_norm;
    MPI_Comm comm;
    HYPRE_Int print_level;
 
@@ -2558,6 +2657,9 @@ void HypreGMRES::Mult(const HypreParVector &b, HypreParVector &x) const
    }
 
    HYPRE_ParCSRGMRESSolve(gmres_solver, *A, b, x);
+   HYPRE_ParCSRGMRESGetNumIterations(gmres_solver, &num_iterations);
+   HYPRE_ParCSRGMRESGetFinalRelativeResidualNorm(gmres_solver,
+                                                 &final_res_norm);
 
    if (print_level > 0)
    {
@@ -2566,10 +2668,6 @@ void HypreGMRES::Mult(const HypreParVector &b, HypreParVector &x) const
       hypre_FinalizeTiming(time_index);
       hypre_ClearTiming();
 
-      HYPRE_ParCSRGMRESGetNumIterations(gmres_solver, &num_iterations);
-      HYPRE_ParCSRGMRESGetFinalRelativeResidualNorm(gmres_solver,
-                                                    &final_res_norm);
-
       MPI_Comm_rank(comm, &myid);
 
       if (myid == 0)
@@ -2781,6 +2879,72 @@ HypreBoomerAMG::HypreBoomerAMG(HypreParMatrix &A) : HypreSolver(&A)
    SetDefaultOptions();
 }
 
+void HypreBoomerAMG::Mult(const HypreParVector &b, HypreParVector &x) const
+{
+   int myid;
+   HYPRE_Int time_index = 0;
+   MPI_Comm comm;
+   HYPRE_Int print_level;
+
+   HYPRE_BoomerAMGGetPrintLevel(amg_precond, &print_level);
+
+   HYPRE_ParCSRMatrixGetComm(*A, &comm);
+
+   if (!setup_called)
+   {
+      if (print_level > 0)
+      {
+         time_index = hypre_InitializeTiming("BoomerAMG Setup");
+         hypre_BeginTiming(time_index);
+      }
+
+      HYPRE_BoomerAMGSetup(amg_precond, *A, b, x);
+      setup_called = 1;
+
+      if (print_level > 0)
+      {
+         hypre_EndTiming(time_index);
+         hypre_PrintTiming("Setup phase times", comm);
+         hypre_FinalizeTiming(time_index);
+         hypre_ClearTiming();
+      }
+   }
+
+   if (print_level > 0)
+   {
+      time_index = hypre_InitializeTiming("BoomerAMG Solve");
+      hypre_BeginTiming(time_index);
+   }
+
+   if (!iterative_mode)
+   {
+      x = 0.0;
+   }
+
+   HYPRE_BoomerAMGSolve(amg_precond, *A, b, x);
+   HYPRE_BoomerAMGGetNumIterations(amg_precond, &num_iterations);
+   HYPRE_BoomerAMGGetFinalRelativeResidualNorm(amg_precond,
+                                               &final_res_norm);
+
+   if (print_level > 0)
+   {
+      hypre_EndTiming(time_index);
+      hypre_PrintTiming("Solve phase times", comm);
+      hypre_FinalizeTiming(time_index);
+      hypre_ClearTiming();
+
+      MPI_Comm_rank(comm, &myid);
+
+      if (myid == 0)
+      {
+         mfem::out << "BoomerAMG Iterations = " << num_iterations << endl
+                   << "Final Relative Residual Norm = " << final_res_norm
+                   << endl;
+      }
+   }
+}
+
+
 void HypreBoomerAMG::SetDefaultOptions()
 {
    // AMG coarsening options:
@@ -3028,6 +3192,250 @@ void HypreBoomerAMG::SetElasticityOptions(ParFiniteElementSpace *fespace)
    error_mode = IGNORE_HYPRE_ERRORS;
 }
 
+void HypreBoomerAMG::SetCoord(int coord_dim, float *coord)
+{
+   HYPRE_BoomerAMGSetPlotGrids (amg_precond, 1);
+   //HYPRE_BoomerAMGSetPlotFileName (amg_precond, plot_file_name);
+   HYPRE_BoomerAMGSetCoordDim (amg_precond, coord_dim);
+   HYPRE_BoomerAMGSetCoordinates (amg_precond, coord);
+}
+
+
+void HypreBoomerAMG::SetLAIROptions(int distance,
+                                    std::string prerelax,
+                                    std::string postrelax,
+                                    double strength_tolC,
+                                    double strength_tolR,
+                                    double filter_tolR,
+                                    int interp_type,
+                                    int relax_type,
+                                    double filterA_tol,
+                                    int splitting,
+                                    int blksize,
+                                    int Sabs)
+{
+   int ns_down, ns_up, ns_coarse;
+   if (distance > 0)
+   {
+      ns_down = prerelax.length();
+      ns_up = postrelax.length();
+      ns_coarse = 1;
+      std::string F("F");
+      std::string C("C");
+      std::string A("A");
+
+      // Array to store relaxation scheme and pass to Hypre
+      int **grid_relax_points = (int **) malloc(4*sizeof(int *));
+      grid_relax_points[0] = NULL;
+      grid_relax_points[1] = (int *) malloc(sizeof(int)*ns_down);
+      grid_relax_points[2] = (int *) malloc(sizeof(int)*ns_up);
+      grid_relax_points[3] = (int *) malloc(sizeof(int));
+      grid_relax_points[3][0] = 0;
+
+      // set down relax scheme
+      for (unsigned int i = 0; i<ns_down; i++)
+      {
+         if (prerelax.compare(i,1,F) == 0)
+         {
+            grid_relax_points[1][i] = -1;
+         }
+         else if (prerelax.compare(i,1,C) == 0)
+         {
+            grid_relax_points[1][i] = 1;
+         }
+         else if (prerelax.compare(i,1,A) == 0)
+         {
+            grid_relax_points[1][i] = 0;
+         }
+      }
+
+      // set up relax scheme
+      for (unsigned int i = 0; i<ns_up; i++)
+      {
+         if (postrelax.compare(i,1,F) == 0)
+         {
+            grid_relax_points[2][i] = -1;
+         }
+         else if (postrelax.compare(i,1,C) == 0)
+         {
+            grid_relax_points[2][i] = 1;
+         }
+         else if (postrelax.compare(i,1,A) == 0)
+         {
+            grid_relax_points[2][i] = 0;
+         }
+      }
+
+      HYPRE_BoomerAMGSetRestriction(amg_precond, distance);
+
+      HYPRE_BoomerAMGSetGridRelaxPoints(amg_precond, grid_relax_points);
+
+      HYPRE_BoomerAMGSetInterpType(amg_precond, interp_type);
+   }
+
+   //HYPRE_BoomerAMGSetMaxRowSum(amg_precond, 0.8);
+   if (Sabs)
+   {
+      HYPRE_BoomerAMGSetSabs(amg_precond, Sabs);
+   }
+
+   if (blksize > 0)
+   {
+      HYPRE_BoomerAMGSetNumFunctions(amg_precond, blksize);
+      HYPRE_BoomerAMGSetNodal(amg_precond, 1);
+      //HYPRE_BoomerAMGSetNodalLevels(amg_precond, 1);
+   }
+
+   HYPRE_BoomerAMGSetCoarsenType(amg_precond, splitting);
+
+   /* does not support aggressive coarsening */
+   HYPRE_BoomerAMGSetAggNumLevels(amg_precond, 0);
+
+   HYPRE_BoomerAMGSetStrongThreshold(amg_precond, strength_tolC);
+
+   if (distance > 0)
+   {
+      HYPRE_BoomerAMGSetStrongThresholdR(amg_precond, strength_tolR);
+      HYPRE_BoomerAMGSetFilterThresholdR(amg_precond, filter_tolR);
+   }
+
+   if (relax_type > -1)
+   {
+      HYPRE_BoomerAMGSetRelaxType(amg_precond, relax_type);
+   }
+
+   if (distance > 0)
+   {
+      HYPRE_BoomerAMGSetCycleNumSweeps(amg_precond, ns_coarse, 3);
+      HYPRE_BoomerAMGSetCycleNumSweeps(amg_precond, ns_down,   1);
+      HYPRE_BoomerAMGSetCycleNumSweeps(amg_precond, ns_up,     2);
+
+      HYPRE_BoomerAMGSetADropTol(amg_precond, filterA_tol);
+      /* type = -1: drop based on row inf-norm */
+      HYPRE_BoomerAMGSetADropType(amg_precond, -1);
+   }
+
+   //HYPRE_BoomerAMGSetMaxCoarseSize(amg_precond, 1000);
+}
+
+
+void HypreBoomerAMG::SetNAIROptions(int neumann_degree,
+                                    std::string prerelax,
+                                    std::string postrelax,
+                                    double strength_tolC,
+                                    double strength_tolR,
+                                    double filter_tolR,
+                                    int interp_type,
+                                    int relax_type,
+                                    double filterA_tol,
+                                    int splitting,
+                                    int blksize,
+                                    int Sabs)
+{
+   int ns_down, ns_up, ns_coarse;
+   if (neumann_degree > 0)
+   {
+      ns_down = prerelax.length();
+      ns_up = postrelax.length();
+      ns_coarse = 1;
+      std::string F("F");
+      std::string C("C");
+      std::string A("A");
+
+      // Array to store relaxation scheme and pass to Hypre
+      int **grid_relax_points = (int **) malloc(4*sizeof(int *));
+      grid_relax_points[0] = NULL;
+      grid_relax_points[1] = (int *) malloc(sizeof(int)*ns_down);
+      grid_relax_points[2] = (int *) malloc(sizeof(int)*ns_up);
+      grid_relax_points[3] = (int *) malloc(sizeof(int));
+      grid_relax_points[3][0] = 0;
+
+      // set down relax scheme
+      for (unsigned int i = 0; i<ns_down; i++)
+      {
+         if (prerelax.compare(i,1,F) == 0)
+         {
+            grid_relax_points[1][i] = -1;
+         }
+         else if (prerelax.compare(i,1,C) == 0)
+         {
+            grid_relax_points[1][i] = 1;
+         }
+         else if (prerelax.compare(i,1,A) == 0)
+         {
+            grid_relax_points[1][i] = 0;
+         }
+      }
+
+      // set up relax scheme
+      for (unsigned int i = 0; i<ns_up; i++)
+      {
+         if (postrelax.compare(i,1,F) == 0)
+         {
+            grid_relax_points[2][i] = -1;
+         }
+         else if (postrelax.compare(i,1,C) == 0)
+         {
+            grid_relax_points[2][i] = 1;
+         }
+         else if (postrelax.compare(i,1,A) == 0)
+         {
+            grid_relax_points[2][i] = 0;
+         }
+      }
+
+      HYPRE_BoomerAMGSetRestriction(amg_precond, 3+neumann_degree);
+
+      HYPRE_BoomerAMGSetGridRelaxPoints(amg_precond, grid_relax_points);
+
+      HYPRE_BoomerAMGSetInterpType(amg_precond, interp_type);
+   }
+
+   //HYPRE_BoomerAMGSetMaxRowSum(amg_precond, 0.8);
+   if (Sabs)
+   {
+      HYPRE_BoomerAMGSetSabs(amg_precond, Sabs);
+   }
+
+   if (blksize > 0)
+   {
+      HYPRE_BoomerAMGSetNumFunctions(amg_precond, blksize);
+      HYPRE_BoomerAMGSetNodal(amg_precond, 1);
+      //HYPRE_BoomerAMGSetNodalLevels(amg_precond, 1);
+   }
+
+   HYPRE_BoomerAMGSetCoarsenType(amg_precond, splitting);
+
+   /* does not support aggressive coarsening */
+   HYPRE_BoomerAMGSetAggNumLevels(amg_precond, 0);
+
+   HYPRE_BoomerAMGSetStrongThreshold(amg_precond, strength_tolC);
+
+   if (neumann_degree > 0)
+   {
+      HYPRE_BoomerAMGSetStrongThresholdR(amg_precond, strength_tolR);
+   }
+
+   if (relax_type > -1)
+   {
+      HYPRE_BoomerAMGSetRelaxType(amg_precond, relax_type);
+   }
+
+   if (neumann_degree > 0)
+   {
+      HYPRE_BoomerAMGSetCycleNumSweeps(amg_precond, ns_coarse, 3);
+      HYPRE_BoomerAMGSetCycleNumSweeps(amg_precond, ns_down,   1);
+      HYPRE_BoomerAMGSetCycleNumSweeps(amg_precond, ns_up,     2);
+
+      HYPRE_BoomerAMGSetADropTol(amg_precond, filterA_tol);
+      /* type = -1: drop based on row inf-norm */
+      HYPRE_BoomerAMGSetADropType(amg_precond, -1);
+   }
+
+   //HYPRE_BoomerAMGSetMaxCoarseSize(amg_precond, 1000);
+}
+
+
 HypreBoomerAMG::~HypreBoomerAMG()
 {
    for (int i = 0; i < rbms.Size(); i++)
@@ -4091,4 +4499,4 @@ HypreAME::StealEigenvectors()
 
 }
 
-#endif
+#endif

linalg/hypre.hpp

@@ -84,6 +84,9 @@ private:
    inline void _SetDataAndSize_();
 
 public:
+
+   HypreParVector() {}
+
    /** @brief Creates vector with given global size and parallel partitioning of
        the rows/columns given by @a col. */
    /** @anchor hypre_partitioning_descr
@@ -116,9 +119,9 @@ public:
    /// MPI communicator
    MPI_Comm GetComm() { return x->comm; }
 
-   /// Returns the parallel row/column partitioning
-   /** See @ref hypre_partitioning_descr "here" for a description of the
-       partitioning array. */
+   void WrapHypreParVector(hypre_ParVector *y);
+
+   /// Returns the row partitioning
    inline HYPRE_Int *Partitioning() { return x->partitioning; }
 
    /// Returns the global number of rows
@@ -234,22 +237,25 @@ public:
    /// An empty matrix to be used as a reference to an existing matrix
    HypreParMatrix();
 
-   /// Converts hypre's format to HypreParMatrix
-   /** If @a owner is false, ownership of @a a is not transferred */
-   explicit HypreParMatrix(hypre_ParCSRMatrix *a, bool owner = true)
+   void WrapHypreParCSRMatrix(hypre_ParCSRMatrix *a, bool owner = true)
    {
-      Init();
       A = a;
       if (!owner) { ParCSROwner = 0; }
       height = GetNumRows();
       width = GetNumCols();
    }
 
-   /// Creates block-diagonal square parallel matrix.
-   /** Diagonal is given by @a diag which must be in CSR format (finalized). The
-       new HypreParMatrix does not take ownership of any of the input arrays.
-       See @ref hypre_partitioning_descr "here" for a description of the row
-       partitioning array @a row_starts.
+   /// Converts hypre's format to HypreParMatrix
+   /** If @a owner is false, ownership of @a a is not transferred */
+   explicit HypreParMatrix(hypre_ParCSRMatrix *a, bool owner = true)
+   {
+      Init();
+      WrapHypreParCSRMatrix(a, owner);
+   }
+
+   /** Creates block-diagonal square parallel matrix. Diagonal is given by diag
+       which must be in CSR format (finalized). The new HypreParMatrix does not
+       take ownership of any of the input arrays.
 
        @warning The ordering of the columns in each row in @a *diag may be
        changed by this constructor to ensure that the first entry in each row is
@@ -334,6 +340,7 @@ public:
 
    /// MPI communicator
    MPI_Comm GetComm() const { return A->comm; }
+   void BuildComm() const { hypre_MatvecCommPkgCreate(A); }
 
    /// Typecasting to hypre's hypre_ParCSRMatrix*
    operator hypre_ParCSRMatrix*() const { return A; }
@@ -393,6 +400,8 @@ public:
    void GetDiag(SparseMatrix &diag) const;
    /// Get the local off-diagonal block. NOTE: 'offd' will not own any data.
    void GetOffd(SparseMatrix &offd, HYPRE_Int* &cmap) const;
+   /// Get on-processor rows as CSR matrix.
+   void GetProcRows(SparseMatrix &colCSRMat);
 
    /** Split the matrix into M x N equally sized blocks of parallel matrices.
        The size of 'blocks' must already be set to M x N. */
@@ -403,6 +412,11 @@ public:
    /// Returns the transpose of *this
    HypreParMatrix * Transpose() const;
 
+   /** Returns principle submatrix given by array of indices of connections
+       with relative size > \@ threshold in *this. */
+   HypreParMatrix * ExtractSubmatrix(Array<int> &indices,
+                                     double threshhold=0.0) const;
+
    /// Returns the number of rows in the diagonal block of the ParCSRMatrix
    int GetNumRows() const
    {
@@ -549,16 +563,22 @@ public:
    Type GetType() const { return Hypre_ParCSR; }
 };
 
+int BlockInvScal(const HypreParMatrix *A, HypreParMatrix *C,
+                 const Vector *b, HypreParVector *d, int block, int job);
+
 /** @brief Return a new matrix `C = alpha*A + beta*B`, assuming that both `A`
     and `B` use the same row and column partitions and the same `col_map_offd`
     arrays. */
 HypreParMatrix *Add(double alpha, const HypreParMatrix &A,
                     double beta,  const HypreParMatrix &B);
+HypreParMatrix *HypreParMatrixAdd(double alpha, const HypreParMatrix &A,
+                                  double beta,  const HypreParMatrix &B);
 
 /** Returns the matrix @a A * @a B. Returned matrix does not necessarily own
     row or column starts unless the bool @a own_matrix is set to true. */
 HypreParMatrix * ParMult(const HypreParMatrix *A, const HypreParMatrix *B,
                          bool own_matrix = false);
+
 /// Returns the matrix A + B
 /** It is assumed that both matrices use the same row and column partitions and
     the same col_map_offd arrays. */
@@ -637,7 +657,7 @@ public:
        1001 = Taubin polynomial smoother
        1002 = FIR polynomial smoother. */
    enum Type { Jacobi = 0, l1Jacobi = 1, l1GS = 2, l1GStr = 4, lumpedJacobi = 5,
-               GS = 6, Chebyshev = 16, Taubin = 1001, FIR = 1002
+               GS = 6, TS = 10, Chebyshev = 16, Taubin = 1001, FIR = 1002
              };
 
    HypreSmoother();
@@ -648,7 +668,7 @@ public:
                  double poly_fraction = .3);
 
    /// Set the relaxation type and number of sweeps
-   void SetType(HypreSmoother::Type type, int relax_times = 1);
+   void SetType(int type, int relax_times = 1);
    /// Set SOR-related parameters
    void SetSOROptions(double relax_weight, double omega);
    /// Set parameters for polynomial smoothing
@@ -701,6 +721,8 @@ protected:
 
    /// Was hypre's Setup function called already?
    mutable int setup_called;
+   mutable HYPRE_Int num_iterations;
+   mutable double final_res_norm;
 
    /// How to treat hypre errors.
    mutable ErrorMode error_mode;
@@ -710,6 +732,9 @@ public:
 
    HypreSolver(HypreParMatrix *_A);
 
+   int GetNumIterations() const { return num_iterations; }
+   double GetFinalNorm() const { return final_res_norm; }
+
    /// Typecast to HYPRE_Solver -- return the solver
    virtual operator HYPRE_Solver() const = 0;
 
@@ -739,6 +764,25 @@ public:
    virtual ~HypreSolver();
 };
 
+
+/// Abstract class for hypre's solvers and preconditioners
+class HypreTriSolve : public HypreSolver
+{
+public:
+   HypreTriSolve() : HypreSolver() { }
+   explicit HypreTriSolve(HypreParMatrix &A) : HypreSolver(&A) { }
+   virtual operator HYPRE_Solver() const { return NULL; }
+
+   virtual HYPRE_PtrToParSolverFcn SetupFcn() const
+   { return (HYPRE_PtrToParSolverFcn) HYPRE_ParCSROnProcTriSetup; }
+   virtual HYPRE_PtrToParSolverFcn SolveFcn() const
+   { return (HYPRE_PtrToParSolverFcn) HYPRE_ParCSROnProcTriSolve; }
+
+   HypreParMatrix* GetData() { return A; }
+   virtual ~HypreTriSolve() { }
+};
+
+
 /// PCG solver in hypre
 class HyprePCG : public HypreSolver
 {
@@ -813,6 +857,7 @@ public:
    virtual void SetOperator(const Operator &op);
 
    void SetTol(double tol);
+   void SetAbsTol(double tol);
    void SetMaxIter(int max_iter);
    void SetKDim(int dim);
    void SetLogging(int logging);
@@ -989,9 +1034,83 @@ public:
        As with SetSystemsOptions(), this solver assumes Ordering::byVDIM. */
    void SetElasticityOptions(ParFiniteElementSpace *fespace);
 
+   /* distance parameter takes on values {1,2,15} for lAIR, meaning R is built using
+   distance 1 neighbors, distance two neighbors, or distance two on processor and
+   distance 1 off processor (i.e., distance 1.5 --> 15).        */
+   void SetLAIROptions(int distance=15,  std::string prerelax="",
+                       std::string postrelax="FFC", double strength_tol=0.1,
+                       double strength_tolR=0.01, double filter_tolR=0.0,
+                       int interp_type=100, int relax_type=3, double filterA_tol=0.0,
+                       int splitting=6, int blksize=0, int Sabs=0);
+
+   void SetNAIROptions(int neumann_degree=2,  std::string prerelax="A",
+                       std::string postrelax="F", double strength_tol=0.1,
+                       double strength_tolR=0.01, double filter_tolR=0.0,
+                       int interp_type=100, int relax_type=10, double filterA_tol=0.0,
+                       int splitting=6, int blksize=0, int Sabs=0);
+
+   void SetCoord(int dim, float *coord);
+
    void SetPrintLevel(int print_level)
    { HYPRE_BoomerAMGSetPrintLevel(amg_precond, print_level); }
 
+   void SetMaxIter(int max_iter)
+   { HYPRE_BoomerAMGSetMaxIter(amg_precond, max_iter); }
+
+   void SetMaxLevels(int max_levels)
+   { HYPRE_BoomerAMGSetMaxLevels(amg_precond, max_levels); }
+
+   void SetTol(double tol)
+   { HYPRE_BoomerAMGSetTol(amg_precond, tol); }
+
+   void SetStrengthThresh(double strength)
+   { HYPRE_BoomerAMGSetStrongThreshold(amg_precond, strength); }
+
+   void SetStrengthThreshR(double strengthR)
+   { HYPRE_BoomerAMGSetStrongThresholdR(amg_precond, strengthR); }
+
+   void SetFilterThreshR(double filterR)
+   { HYPRE_BoomerAMGSetFilterThresholdR(amg_precond, filterR); }
+
+   void SetInterpolation(int interp_type)
+   { HYPRE_BoomerAMGSetInterpType(amg_precond, interp_type); }
+
+   void SetRestriction(int restrict_type)
+   { HYPRE_BoomerAMGSetRestriction(amg_precond, restrict_type); }
+
+   void SetCoarsening(int coarsen_type)
+   { HYPRE_BoomerAMGSetCoarsenType(amg_precond, coarsen_type); }
+
+   void SetRelaxType(int relax_type)
+   { HYPRE_BoomerAMGSetRelaxType(amg_precond, relax_type); }
+
+   void SetRelaxCycle(int prerelax, int postrelax)
+   {
+      HYPRE_BoomerAMGSetCycleNumSweeps(amg_precond, prerelax,  1);
+      HYPRE_BoomerAMGSetCycleNumSweeps(amg_precond, postrelax, 2);
+   }
+
+   void GetNumIterations(int &num_it)
+   { HYPRE_BoomerAMGGetNumIterations(amg_precond, &num_it); }
+
+   void SetCycleType(int cycle_type)
+   { HYPRE_BoomerAMGSetCycleType(amg_precond, cycle_type); }
+
+   void SetNodal(int blocksize)
+   {
+      HYPRE_BoomerAMGSetNumFunctions(amg_precond, blocksize);
+      HYPRE_BoomerAMGSetNodal(amg_precond, 1);
+   }
+
+   void SetAggressiveCoarsening(int num_levels)
+   { HYPRE_BoomerAMGSetAggNumLevels(amg_precond, num_levels); }
+
+   void SetTriangular()
+   { HYPRE_BoomerAMGSetIsTriangular(amg_precond, 1); }
+
+   void SetGMRESSwitchR(int gmres_switch)
+   { HYPRE_BoomerAMGSetGMRESSwitchR(amg_precond, gmres_switch); }
+
    /// The typecast to HYPRE_Solver returns the internal amg_precond
    virtual operator HYPRE_Solver() const { return amg_precond; }
 
@@ -1000,6 +1119,9 @@ public:
    virtual HYPRE_PtrToParSolverFcn SolveFcn() const
    { return (HYPRE_PtrToParSolverFcn) HYPRE_BoomerAMGSolve; }
 
+   virtual void Mult (const HypreParVector &b, HypreParVector &x) const;
+   using HypreSolver::Mult;
+
    virtual ~HypreBoomerAMG();
 };
 

linalg/solvers.cpp

@@ -232,23 +232,23 @@ void SLISolver::Mult(const Vector &b, Vector &x) const
    if (prec)
    {
       prec->Mult(r, z); // z = B r
-      nom0 = nom = Dot(z, r);
+      nom0 = nom = sqrt(Dot(z, z));
    }
    else
    {
-      nom0 = nom = Dot(r, r);
+      nom0 = nom = sqrt(Dot(r, r));
    }
 
    if (print_level == 1)
-      mfem::out << "   Iteration : " << setw(3) << 0 << "  (B r, r) = "
+      mfem::out << "   Iteration : " << setw(3) << 0 << "  ||Br|| = "
                 << nom << '\n';
 
-   r0 = std::max(nom*rel_tol*rel_tol, abs_tol*abs_tol);
+   r0 = std::max(nom*rel_tol, abs_tol);
    if (nom <= r0)
    {
       converged = 1;
       final_iter = 0;
-      final_norm = sqrt(nom);
+      final_norm = nom;
       return;
    }
 
@@ -272,16 +272,16 @@ void SLISolver::Mult(const Vector &b, Vector &x) const
       if (prec)
       {
          prec->Mult(r, z); //  z = B r
-         nom = Dot(z, r);
+         nom = sqrt(Dot(z, z));
       }
       else
       {
-         nom = Dot(r, r);
+         nom = sqrt(Dot(r, r));
       }
 
-      cf = sqrt(nom/nomold);
+      cf = nom/nomold;
       if (print_level == 1)
-         mfem::out << "   Iteration : " << setw(3) << i << "  (B r, r) = "
+         mfem::out << "   Iteration : " << setw(3) << i << "  ||Br|| = "
                    << nom << "\tConv. rate: " << cf << '\n';
       nomold = nom;
 
@@ -291,8 +291,8 @@ void SLISolver::Mult(const Vector &b, Vector &x) const
             mfem::out << "Number of SLI iterations: " << i << '\n'
                       << "Conv. rate: " << cf << '\n';
          else if (print_level == 3)
-            mfem::out << "(B r_0, r_0) = " << nom0 << '\n'
-                      << "(B r_N, r_N) = " << nom << '\n'
+            mfem::out << "||Br_0|| = " << nom0 << '\n'
+                      << "||Br_N|| = " << nom << '\n'
                       << "Number of SLI iterations: " << i << '\n';
          converged = 1;
          final_iter = i;
@@ -308,16 +308,16 @@ void SLISolver::Mult(const Vector &b, Vector &x) const
    if (print_level >= 0 && !converged)
    {
       mfem::err << "SLI: No convergence!" << '\n';
-      mfem::out << "(B r_0, r_0) = " << nom0 << '\n'
-                << "(B r_N, r_N) = " << nom << '\n'
+      mfem::out << "||Br_0|| = " << nom0 << '\n'
+                << "||Br_N|| = " << nom << '\n'
                 << "Number of SLI iterations: " << final_iter << '\n';
    }
    if (print_level >= 1 || (print_level >= 0 && !converged))
    {
       mfem::out << "Average reduction factor = "
-                << pow (nom/nom0, 0.5/final_iter) << '\n';
+                << pow (nom/nom0, 1.0/final_iter) << '\n';
    }
-   final_norm = sqrt(nom);
+   final_norm = nom;
 }
 
 void SLI(const Operator &A, const Vector &b, Vector &x,

linalg/vector.cpp

@@ -495,6 +495,22 @@ void Vector::GetSubVector(const Array<int> &dofs, double *elem_data) const
    }
 }
 
+// ADDED //
+void Vector::GetSubVector(int index_low, int index_high, Vector &elemvect) const
+{
+    int i, j, n = index_high - index_low;
+    
+    elemvect.SetSize (n);
+    
+    int k = 0;
+    for (i = index_low; i < index_high; i++)
+    {
+        elemvect(k) = data[i];
+        k += 1;
+    }
+}
+// ADDED //
+
 void Vector::SetSubVector(const Array<int> &dofs, const double value)
 {
    const bool use_dev = dofs.UseDevice();
@@ -561,6 +577,36 @@ void Vector::SetSubVector(const Array<int> &dofs, double *elem_data)
    }
 }
 
+
+
+// ADDED //
+void Vector::SetSubVector(int index_low, int index_high, Vector &elemvect)
+{
+    int i, j, n = index_high - index_low;
+    
+    int k = 0;
+    for (i = index_low; i < index_high; i++)
+    {
+        data[i] = elemvect(k);
+        k += 1;
+    }
+}
+
+
+void Vector::AddElementVector(int index_low, int index_high, double c, Vector &elemvect)
+{
+    int i, j;
+    
+    int k = 0;
+    for (i = index_low; i < index_high; i++)
+    {
+        data[i] += c * elemvect(k);
+        k += 1;
+    }
+}
+// ADDED //
+
+
 void Vector::AddElementVector(const Array<int> &dofs, const Vector &elemvect)
 {
    MFEM_ASSERT(dofs.Size() == elemvect.Size(), "Size mismatch: "

linalg/vector.hpp

@@ -281,17 +281,29 @@ public:
    void GetSubVector(const Array<int> &dofs, Vector &elemvect) const;
    void GetSubVector(const Array<int> &dofs, double *elem_data) const;
 
+   // ADDED //
+   void GetSubVector(int index_lo, int index_high, Vector &elemvect) const;
+   // ADDED //
+
    /// Set the entries listed in `dofs` to the given `value`.
    void SetSubVector(const Array<int> &dofs, const double value);
    void SetSubVector(const Array<int> &dofs, const Vector &elemvect);
    void SetSubVector(const Array<int> &dofs, double *elem_data);
 
+   // ADDED //
+   void SetSubVector(int index_low, int index_high, Vector &elemvect);
+   // ADDED //
+
    /// Add (element) subvector to the vector.
    void AddElementVector(const Array<int> & dofs, const Vector & elemvect);
    void AddElementVector(const Array<int> & dofs, double *elem_data);
    void AddElementVector(const Array<int> & dofs, const double a,
                          const Vector & elemvect);
 
+   // ADDED //
+   void AddElementVector(int index_low, int index_high, double c, Vector &elemvect);
+   // ADDED //
+
    /// Set all vector entries NOT in the 'dofs' array to the given 'val'.
    void SetSubVectorComplement(const Array<int> &dofs, const double val);
 

mesh/mesh.cpp

@@ -926,6 +926,82 @@ FaceElementTransformations *Mesh::GetFaceElementTransformations(int FaceNo,
    return &FaceElemTr;
 }
 
+
+
+// ADDED //
+void Mesh::GetFaceElementTransformations(int FaceNo, 
+                                         FaceElementTransformations &FaceElemTr,
+                                         IsoparametricTransformation &Transformation, 
+                                         IsoparametricTransformation &Transformation2,
+                                         IsoparametricTransformation &FTr)
+{
+
+   int mask = 31;
+   FaceInfo &face_info = faces_info[FaceNo];
+
+   FaceElemTr.Elem1 = NULL;
+   FaceElemTr.Elem2 = NULL;
+
+   // setup the transformation for the first element
+   FaceElemTr.Elem1No = face_info.Elem1No;
+   if (mask & 1)
+   {
+      GetElementTransformation(FaceElemTr.Elem1No, &Transformation);
+      FaceElemTr.Elem1 = &Transformation;
+   }
+
+   //  setup the transformation for the second element
+   //     return NULL in the Elem2 field if there's no second element, i.e.
+   //     the face is on the "boundary"
+   FaceElemTr.Elem2No = face_info.Elem2No;
+   if ((mask & 2) && FaceElemTr.Elem2No >= 0)
+   {
+#ifdef MFEM_DEBUG
+      if (NURBSext && (mask & 1)) { MFEM_ABORT("NURBS mesh not supported!"); }
+#endif
+      GetElementTransformation(FaceElemTr.Elem2No, &Transformation2);
+      FaceElemTr.Elem2 = &Transformation2;
+   }
+
+   // setup the face transformation
+   FaceElemTr.FaceGeom = GetFaceGeometryType(FaceNo);
+   GetFaceTransformation(FaceNo, &FTr);
+   FaceElemTr.Face = &FTr;
+
+   // setup Loc1 & Loc2
+   int face_type = GetFaceElementType(FaceNo);
+   if (mask & 4)
+   {
+      int elem_type = GetElementType(face_info.Elem1No);
+      GetLocalFaceTransformation(face_type, elem_type,
+                                 FaceElemTr.Loc1.Transf, face_info.Elem1Inf);
+   }
+   if ((mask & 8) && FaceElemTr.Elem2No >= 0)
+   {
+      int elem_type = GetElementType(face_info.Elem2No);
+      GetLocalFaceTransformation(face_type, elem_type,
+                                 FaceElemTr.Loc2.Transf, face_info.Elem2Inf);
+
+      // NC meshes: prepend slave edge/face transformation to Loc2
+      if (Nonconforming() && IsSlaveFace(face_info))
+      {
+         ApplyLocalSlaveTransformation(FaceElemTr.Loc2.Transf, face_info);
+
+         if (face_type == Element::SEGMENT)
+         {
+            // flip Loc2 to match Loc1 and Face
+            DenseMatrix &pm = FaceElemTr.Loc2.Transf.GetPointMat();
+            std::swap(pm(0,0), pm(0,1));
+            std::swap(pm(1,0), pm(1,1));
+         }
+      }
+   }
+
+}
+// ADDED //
+
+
+
 bool Mesh::IsSlaveFace(const FaceInfo &fi) const
 {
    return fi.NCFace >= 0 && nc_faces_info[fi.NCFace].Slave;
@@ -971,6 +1047,45 @@ FaceElementTransformations *Mesh::GetBdrFaceTransformations(int BdrElemNo)
    return tr;
 }
 
+
+// ADDED //
+bool Mesh::GetBdrFaceTransformations(int BdrElemNo,
+                                     FaceElementTransformations &tr,
+                                     IsoparametricTransformation &Transformation, 
+                                     IsoparametricTransformation &Transformation2,
+                                     IsoparametricTransformation &FTr)
+{
+   int fn;
+   bool is_not_null = true;
+   if (Dim == 3)
+   {
+      fn = be_to_face[BdrElemNo];
+   }
+   else if (Dim == 2)
+   {
+      fn = be_to_edge[BdrElemNo];
+   }
+   else
+   {
+      fn = boundary[BdrElemNo]->GetVertices()[0];
+   }
+   // Check if the face is interior, shared, or non-conforming.
+   if (FaceIsTrueInterior(fn) || faces_info[fn].NCFace >= 0)
+   {
+      is_not_null = false;
+      return is_not_null;
+   }
+   GetFaceElementTransformations(fn, tr, Transformation, 
+                                 Transformation2, FTr);
+   tr.Face->Attribute = boundary[BdrElemNo]->GetAttribute();
+
+   return is_not_null;
+
+}
+// ADDED //
+
+
+
 void Mesh::GetFaceElements(int Face, int *Elem1, int *Elem2) const
 {
    *Elem1 = faces_info[Face].Elem1No;

mesh/mesh.hpp

@@ -971,14 +971,42 @@ public:
    FaceElementTransformations *GetFaceElementTransformations(int FaceNo,
                                                              int mask = 31);
 
+   // ADDED //
+   void GetFaceElementTransformations(int FaceNo, 
+                                      FaceElementTransformations &FaceElemTr,
+                                      IsoparametricTransformation &Transformation, 
+                                      IsoparametricTransformation &Transformation2,
+                                      IsoparametricTransformation &FTr);
+   // ADDED //
+
    FaceElementTransformations *GetInteriorFaceTransformations (int FaceNo)
    {
       if (faces_info[FaceNo].Elem2No < 0) { return NULL; }
       return GetFaceElementTransformations (FaceNo);
    }
 
+   // ADDED //
+   void GetInteriorFaceTransformations (int FaceNo, 
+                                        FaceElementTransformations &FaceElemTr,
+                                        IsoparametricTransformation &Transformation, 
+                                        IsoparametricTransformation &Transformation2,
+                                        IsoparametricTransformation &FTr)
+   {
+      GetFaceElementTransformations (FaceNo, FaceElemTr, 
+                                     Transformation, Transformation2, FTr);
+   }
+   // ADDED //
+
    FaceElementTransformations *GetBdrFaceTransformations (int BdrElemNo);
 
+   // ADDED //
+   bool GetBdrFaceTransformations(int BdrElemNo,
+                                  FaceElementTransformations &tr,
+                                  IsoparametricTransformation &Transformation, 
+                                  IsoparametricTransformation &Transformation2,
+                                  IsoparametricTransformation &FTr);
+   // ADDED //
+
    /// Return true if the given face is interior. @sa FaceIsTrueInterior().
    bool FaceIsInterior(int FaceNo) const
    {

mesh/pmesh.cpp

@@ -314,12 +314,20 @@ int ParMesh::BuildLocalVertices(const mfem::Mesh &mesh,
 
    vertices.SetSize(vert_counter);
 
+
+   // ADDED //
+   vert_local_to_global.SetSize(mesh.GetNV());
+   vert_local_to_global = -1;
+   // ADDED //
    for (int i = 0; i < vert_global_local.Size(); i++)
    {
       if (vert_global_local[i] >= 0)
       {
          vertices[vert_global_local[i]].SetCoords(mesh.SpaceDimension(),
                                                   mesh.GetVertex(i));
+         // ADDED //
+         vert_local_to_global[vert_global_local[i]] = i;
+         // ADDED //
       }
    }
 
@@ -2438,6 +2446,137 @@ GetSharedFaceTransformations(int sf, bool fill2)
    return &FaceElemTr;
 }
 
+
+// ------------------------ ADDED ------------------------ //
+ElementTransformation* ParMesh::GetGhostFaceTransformation(
+                             FaceElementTransformations* FETr, 
+                             IsoparametricTransformation &FaceTransformation,
+                             Element::Type face_type,
+                             Geometry::Type face_geom)
+{
+   // calculate composition of FETr->Loc1 and FETr->Elem1
+   DenseMatrix &face_pm = FaceTransformation.GetPointMat();
+   if (Nodes == NULL)
+   {
+      FETr->Elem1->Transform(FETr->Loc1.Transf.GetPointMat(), face_pm);
+      FaceTransformation.SetFE(GetTransformationFEforElementType(face_type));
+   }
+   else
+   {
+      const FiniteElement* face_el =
+         Nodes->FESpace()->GetTraceElement(FETr->Elem1No, face_geom);
+
+#if 0 // TODO: handle the case of non-interpolatory Nodes
+      DenseMatrix I;
+      face_el->Project(Transformation.GetFE(), FETr->Loc1.Transf, I);
+      MultABt(Transformation.GetPointMat(), I, pm_face);
+#else
+      IntegrationRule eir(face_el->GetDof());
+      FETr->Loc1.Transform(face_el->GetNodes(), eir);
+      Nodes->GetVectorValues(*FETr->Elem1, eir, face_pm);
+#endif
+      FaceTransformation.SetFE(face_el);
+   }
+   FaceTransformation.FinalizeTransformation();
+   return &FaceTransformation;
+}
+
+FaceElementTransformations *ParMesh::
+GetSharedFaceTransformations(int sf, 
+                             FaceElementTransformations &FaceElemTr,
+                             IsoparametricTransformation &Transformation, 
+                             IsoparametricTransformation &Transformation2,
+                             IsoparametricTransformation &FTr)
+{
+
+   bool fill2 = true;
+
+   int FaceNo = GetSharedFace(sf);
+
+   FaceInfo &face_info = faces_info[FaceNo];
+
+   bool is_slave = Nonconforming() && IsSlaveFace(face_info);
+   bool is_ghost = Nonconforming() && FaceNo >= GetNumFaces();
+
+   NCFaceInfo* nc_info = NULL;
+   if (is_slave) { nc_info = &nc_faces_info[face_info.NCFace]; }
+
+   int local_face = is_ghost ? nc_info->MasterFace : FaceNo;
+   Element::Type  face_type = GetFaceElementType(local_face);
+   Geometry::Type face_geom = GetFaceGeometryType(local_face);
+
+   // setup the transformation for the first element
+   FaceElemTr.Elem1No = face_info.Elem1No;
+   GetElementTransformation(FaceElemTr.Elem1No, &Transformation);
+   FaceElemTr.Elem1 = &Transformation;
+
+   // setup the transformation for the second (neighbor) element
+   if (fill2)
+   {
+      FaceElemTr.Elem2No = -1 - face_info.Elem2No;
+      GetFaceNbrElementTransformation(FaceElemTr.Elem2No, &Transformation2);
+      FaceElemTr.Elem2 = &Transformation2;
+   }
+   else
+   {
+      FaceElemTr.Elem2No = -1;
+   }
+
+   // setup the face transformation if the face is not a ghost
+   FaceElemTr.FaceGeom = face_geom;
+   if (!is_ghost)
+   {
+      FaceElemTr.Face = GetFaceTransformation(FaceNo);
+      // NOTE: The above call overwrites FaceElemTr.Loc1
+   }
+
+   // setup Loc1 & Loc2
+   int elem_type = GetElementType(face_info.Elem1No);
+   GetLocalFaceTransformation(face_type, elem_type, FaceElemTr.Loc1.Transf,
+                              face_info.Elem1Inf);
+
+   if (fill2)
+   {
+      elem_type = face_nbr_elements[FaceElemTr.Elem2No]->GetType();
+      GetLocalFaceTransformation(face_type, elem_type, FaceElemTr.Loc2.Transf,
+                                 face_info.Elem2Inf);
+   }
+
+   // adjust Loc1 or Loc2 of the master face if this is a slave face
+   if (is_slave)
+   {
+      // is a ghost slave? -> master not a ghost -> choose Elem1 local transf
+      // not a ghost slave? -> master is a ghost -> choose Elem2 local transf
+      IsoparametricTransformation &loctr =
+         is_ghost ? FaceElemTr.Loc1.Transf : FaceElemTr.Loc2.Transf;
+
+      if (is_ghost || fill2)
+      {
+         ApplyLocalSlaveTransformation(loctr, face_info);
+      }
+
+      if (face_type == Element::SEGMENT && fill2)
+      {
+         // fix slave orientation in 2D: flip Loc2 to match Loc1 and Face
+         DenseMatrix &pm = FaceElemTr.Loc2.Transf.GetPointMat();
+         std::swap(pm(0,0), pm(0,1));
+         std::swap(pm(1,0), pm(1,1));
+      }
+   }
+
+   // for ghost faces we need a special version of GetFaceTransformation
+   if (is_ghost)
+   {
+      FaceElemTr.Face =
+         GetGhostFaceTransformation(&FaceElemTr, face_type, face_geom);
+   }
+
+   return &FaceElemTr;
+}
+// ------------------------ ADDED ------------------------ //
+
+
+
 int ParMesh::GetNSharedFaces() const
 {
    if (Conforming())

mesh/pmesh.hpp

@@ -100,13 +100,16 @@ protected:
    bool DecodeFaceSplittings(HashTable<Hashed2> &v_to_v, const int *v,
                              const Array<unsigned> &codes, int &pos);
 
-   void GetFaceNbrElementTransformation(
-      int i, IsoparametricTransformation *ElTr);
-
    ElementTransformation* GetGhostFaceTransformation(
       FaceElementTransformations* FETr, Element::Type face_type,
       Geometry::Type face_geom);
 
+   // ADDED //
+   ElementTransformation* GetGhostFaceTransformation(FaceElementTransformations* FETr, 
+                                   IsoparametricTransformation &FaceTransformation,
+                                   Element::Type face_type, Geometry::Type face_geom);
+   // ADDED //
+
    /// Update the groups after triangle refinement
    void RefineGroups(const DSTable &v_to_v, int *middle);
 
@@ -244,10 +247,50 @@ public:
    Table            send_face_nbr_elements;
    Table            send_face_nbr_vertices;
 
+   // ADDED //
+   Array<int> shared_face_to_global_face;
+   Array<int> shared_face_to_MPI_rank;
+   Array<int> vert_local_to_global;
+   int elem_local_to_global;
+   Table group_sface; // in 3D, union of group_stria and group_squad
+   // ADDED //
+
    ParNCMesh* pncmesh;
 
    int GetNGroups() const { return gtopo.NGroups(); }
 
+   // ADDED //
+
+   Table const *GetSharedFacesInGroups() 
+   { 
+
+       // determine whether faces are quads or triangles
+       // NOTE: this assumes all mesh elements have the same geometry type
+       Array<int> verts;
+       GetFaceVertices(0, verts);
+       int nv = verts.Size();
+
+       if (Dim == 3 && nv == 3) 
+       {
+           return &group_stria; 
+       } 
+       else if (Dim == 3 && nv == 4) 
+       {
+           return &group_squad; 
+       } 
+       else
+       {
+           return &group_sedge; 
+       }
+   }
+
+   // ADDED //
+
+   // ADDED (Moved from protected) //
+   void GetFaceNbrElementTransformation(
+      int i, IsoparametricTransformation *ElTr);
+   // ADDED //
+
    ///@{ @name These methods require group > 0
    int GroupNVertices(int group) { return group_svert.RowSize(group-1); }
    int GroupNEdges(int group)    { return group_sedge.RowSize(group-1); }
@@ -284,6 +327,14 @@ public:
    FaceElementTransformations *
    GetSharedFaceTransformations(int sf, bool fill2 = true);
 
+   // ADDED //
+   FaceElementTransformations * GetSharedFaceTransformations(int sf, 
+                                     FaceElementTransformations &FaceElemTr,
+                                     IsoparametricTransformation &Transformation, 
+                                     IsoparametricTransformation &Transformation2,
+                                     IsoparametricTransformation &FTr);
+   // ADDED //
+
    /// Return the number of shared faces (3D), edges (2D), vertices (1D)
    int GetNSharedFaces() const;