Squashing commits on ho-trt

examples/ex9p.cpp

@@ -64,6 +64,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:
@@ -71,7 +138,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)
@@ -79,20 +147,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)
@@ -121,10 +209,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
@@ -143,6 +233,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);
@@ -178,7 +270,11 @@ int main(int argc, char *argv[])
    bool adios2 = 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);
 
@@ -212,6 +308,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.");
@@ -307,7 +407,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();
@@ -316,6 +416,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.
@@ -476,11 +579,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; )
@@ -563,6 +674,7 @@ int main(int argc, char *argv[])
    }
 #endif
    delete dc;
+   delete adv;
 
    MPI_Finalize();
    return 0;
@@ -614,6 +726,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

@@ -268,6 +268,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)
@@ -1130,9 +1156,12 @@ MixedBilinearForm::MixedBilinearForm (FiniteElementSpace *tr_fes,
    bbfi = mbf->bbfi;
    tfbfi = mbf->tfbfi;
    btfbfi = mbf->btfbfi;
+   fbfi = mbf->fbfi; 
+   bfbfi = mbf->bfbfi; 
 
    bbfi_marker = mbf->bbfi_marker;
    btfbfi_marker = mbf->btfbfi_marker;
+   bfbfi_marker = mbf->bfbfi_marker; 
 
    assembly = AssemblyLevel::LEGACYFULL;
    ext = NULL;
@@ -1286,6 +1315,24 @@ void MixedBilinearForm::AddBdrTraceFaceIntegrator(BilinearFormIntegrator *bfi,
    btfbfi_marker.Append(&bdr_marker);
 }
 
+void MixedBilinearForm::AddInteriorFaceIntegrator(BilinearFormIntegrator *bfi) 
+{
+	fbfi.Append(bfi); 
+}
+
+void MixedBilinearForm::AddBdrFaceIntegrator(BilinearFormIntegrator *bfi)
+{
+	bfbfi.Append(bfi); 
+	bfbfi_marker.Append(NULL); 
+}
+
+void MixedBilinearForm::AddBdrFaceIntegrator(BilinearFormIntegrator *bfi,
+											 Array<int> &bdr_marker)
+{
+	bfbfi.Append(bfi); 
+	bfbfi_marker.Append(&bdr_marker); 
+}
+
 void MixedBilinearForm::Assemble (int skip_zeros)
 {
    if (ext)
@@ -1457,6 +1504,97 @@ void MixedBilinearForm::Assemble (int skip_zeros)
          }
       }
    }
+
+   if (fbfi.Size()) 
+   {
+      Array<int> te_vdofs2, tr_vdofs2; 
+      FaceElementTransformations *ftr; 
+      const FiniteElement *test_fe1, *test_fe2, *trial_fe1, *trial_fe2; 
+      int nfaces = mesh->GetNumFaces(); 
+      for (int f=0; f<nfaces; f++) 
+      {
+         ftr = mesh->GetInteriorFaceTransformations(f); 
+         if (ftr) 
+         {
+            trial_fes->GetElementVDofs(ftr->Elem1No, tr_vdofs); 
+            trial_fes->GetElementVDofs(ftr->Elem2No, tr_vdofs2); 
+            tr_vdofs.Append(tr_vdofs2); 
+            trial_fe1 = trial_fes->GetFE(ftr->Elem1No); 
+            trial_fe2 = trial_fes->GetFE(ftr->Elem2No); 
+
+            test_fes->GetElementVDofs(ftr->Elem1No, te_vdofs); 
+            test_fes->GetElementVDofs(ftr->Elem2No, te_vdofs2); 
+            te_vdofs.Append(te_vdofs2); 
+            test_fe1 = test_fes->GetFE(ftr->Elem1No); 
+            test_fe2 = test_fes->GetFE(ftr->Elem2No); 
+
+            for (int k=0; k<fbfi.Size(); k++) 
+            {
+               fbfi[k] -> AssembleFaceMatrix(*trial_fe1, *trial_fe2,
+                                             *test_fe1, *test_fe2, 
+                                             *ftr, elemmat); 
+               mat->AddSubMatrix(te_vdofs, tr_vdofs, elemmat, skip_zeros); 
+            }
+         }
+
+      }
+   }
+
+   if (bfbfi.Size())
+   {
+      FaceElementTransformations *tr;
+      const FiniteElement *tr_fe1, *tr_fe2, *te_fe1, *te_fe2;
+
+      // Which boundary attributes need to be processed?
+      Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
+                                 mesh->bdr_attributes.Max() : 0);
+      bdr_attr_marker = 0;
+      for (int k = 0; k < bfbfi.Size(); k++)
+      {
+         if (bfbfi_marker[k] == NULL)
+         {
+            bdr_attr_marker = 1;
+            break;
+         }
+         Array<int> &bdr_marker = *bfbfi_marker[k];
+         MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
+                     "invalid boundary marker for boundary face integrator #"
+                     << k << ", counting from zero");
+         for (int i = 0; i < bdr_attr_marker.Size(); i++)
+         {
+            bdr_attr_marker[i] |= bdr_marker[i];
+         }
+      }
+
+      for (int i = 0; i < trial_fes -> GetNBE(); i++)
+      {
+         const int bdr_attr = mesh->GetBdrAttribute(i);
+         if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
+
+         tr = mesh -> GetBdrFaceTransformations (i);
+         if (tr != NULL)
+         {
+            trial_fes -> GetElementVDofs (tr -> Elem1No, tr_vdofs);
+            tr_fe1 = trial_fes -> GetFE (tr -> Elem1No);
+            // The fe2 object is really a dummy and not used on the boundaries,
+            // but we can't dereference a NULL pointer, and we don't want to
+            // actually make a fake element.
+            tr_fe2 = tr_fe1;
+
+            test_fes -> GetElementVDofs (tr -> Elem1No, te_vdofs); 
+            te_fe1 = test_fes -> GetFE (tr -> Elem1No); 
+            te_fe2 = te_fe1; 
+            for (int k = 0; k < bfbfi.Size(); k++)
+            {
+               if (bfbfi_marker[k] &&
+                   (*bfbfi_marker[k])[bdr_attr-1] == 0) { continue; }
+
+               bfbfi[k] -> AssembleFaceMatrix (*tr_fe1, *tr_fe2, *te_fe1, *te_fe2, *tr, elemmat);
+               mat -> AddSubMatrix (te_vdofs, tr_vdofs, elemmat, skip_zeros);
+            }
+         }
+      }
+   }
 }
 
 void MixedBilinearForm::AssembleDiagonal_ADAt(const Vector &D,
@@ -1766,6 +1904,8 @@ MixedBilinearForm::~MixedBilinearForm()
       for (i = 0; i < bbfi.Size(); i++) { delete bbfi[i]; }
       for (i = 0; i < tfbfi.Size(); i++) { delete tfbfi[i]; }
       for (i = 0; i < btfbfi.Size(); i++) { delete btfbfi[i]; }
+      for (i = 0; i < fbfi.Size(); i++) { delete fbfi[i]; }
+      for (i = 0; i < bfbfi.Size(); i++) { delete bfbfi[i]; }
    }
    delete ext;
 }

fem/bilinearform.hpp

@@ -349,6 +349,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
+
    /// Sets all sparse values of \f$ M \f$ and \f$ M_e \f$ to 'a'.
    void operator=(const double a)
    {
@@ -642,6 +656,13 @@ protected:
    Array<BilinearFormIntegrator*> btfbfi;
    Array<Array<int>*>             btfbfi_marker;///< Entries are not owned.
 
+   // face integrators for vector DG spaces 
+   Array<BilinearFormIntegrator*> fbfi; 
+
+   // boundary face integrators for vector DG spaces
+   Array<BilinearFormIntegrator*> bfbfi; 
+   Array<Array<int>*> bfbfi_marker; 
+
    DenseMatrix elemmat;
    Array<int>  trial_vdofs, test_vdofs;
 
@@ -733,6 +754,11 @@ public:
    void AddBdrTraceFaceIntegrator (BilinearFormIntegrator * bfi,
                                    Array<int> &bdr_marker);
 
+   void AddInteriorFaceIntegrator(BilinearFormIntegrator *bfi); 
+   void AddBdrFaceIntegrator(BilinearFormIntegrator *bfi); 
+   void AddBdrFaceIntegrator(BilinearFormIntegrator *bfi, 
+                              Array<int> &bdr_marker); 
+
    /// Access all integrators added with AddDomainIntegrator().
    Array<BilinearFormIntegrator*> *GetDBFI() { return &dbfi; }
 

fem/bilininteg.cpp

@@ -153,6 +153,15 @@ void BilinearFormIntegrator::AssembleFaceMatrix(
               " Integrator class.");
 }
 
+void BilinearFormIntegrator::AssembleFaceMatrix(
+   const FiniteElement &trial_fe1, const FiniteElement &trial_fe2, 
+   const FiniteElement &test_fe1, const FiniteElement &test_fe2, 
+   FaceElementTransformations &Trans, DenseMatrix &elmat)
+{
+   MFEM_ABORT("AssembleFaceMatrix (mixed form) is not implemented for this"
+              " Integrator class.");
+}
+
 void BilinearFormIntegrator::AssembleElementVector(
    const FiniteElement &el, ElementTransformation &Tr, const Vector &elfun,
    Vector &elvect)
@@ -202,6 +211,15 @@ void TransposeIntegrator::AssembleFaceMatrix (
    elmat.Transpose (bfi_elmat);
 }
 
+void TransposeIntegrator::AssembleFaceMatrix(
+   const FiniteElement &tr_fe1, const FiniteElement &tr_fe2, 
+   const FiniteElement &te_fe1, const FiniteElement &te_fe2, 
+   FaceElementTransformations &T, DenseMatrix &elmat) 
+{
+   bfi -> AssembleFaceMatrix(tr_fe1, tr_fe2, te_fe1, te_fe2, T, bfi_elmat); 
+   elmat.Transpose(bfi_elmat); 
+}
+
 void LumpedIntegrator::AssembleElementMatrix (
    const FiniteElement &el, ElementTransformation &Trans, DenseMatrix &elmat)
 {
@@ -2603,6 +2621,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/bilininteg.hpp

@@ -155,6 +155,13 @@ public:
                                    FaceElementTransformations &Trans,
                                    DenseMatrix &elmat);
 
+   virtual void AssembleFaceMatrix(const FiniteElement &trial_fe1, 
+   								   const FiniteElement &trial_fe2, 
+   								   const FiniteElement &test_fe1, 
+   								   const FiniteElement &test_fe2, 
+   								   FaceElementTransformations &Trans, 
+   								   DenseMatrix &elmat); 
+
    /// @brief Perform the local action of the BilinearFormIntegrator.
    /// Note that the default implementation in the base class is general but not
    /// efficient.
@@ -272,6 +279,9 @@ public:
                                    const FiniteElement &el2,
                                    FaceElementTransformations &Trans,
                                    DenseMatrix &elmat);
+   virtual void AssembleFaceMatrix(const FiniteElement &tr_fe1, const FiniteElement &tr_fe2, 
+                                   const FiniteElement &te_fe1, const FiniteElement &te_fe2, 
+                                   FaceElementTransformations &Trans, DenseMatrix &elmat); 
 
    using BilinearFormIntegrator::AssemblePA;
 

fem/eltrans.cpp

@@ -552,6 +552,64 @@ 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 /////////////////////////////////
+
+
 void FaceElementTransformations::SetIntPoint(const IntegrationPoint *face_ip)
 {
    IsoparametricTransformation::SetIntPoint(face_ip);

fem/eltrans.hpp

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

fem/fe.cpp

@@ -10732,6 +10732,61 @@ void RT_QuadrilateralElement::CalcDivShape(const IntegrationPoint &ip,
       }
 }
 
+void RT_QuadrilateralElement::CalcDShape(const IntegrationPoint &ip, 
+                                         DenseMatrix &gshape) const 
+{
+   const int pp1 = order;
+
+#ifdef MFEM_THREAD_SAFE
+   Vector shape_cx(pp1 + 1), shape_ox(pp1), shape_cy(pp1 + 1), shape_oy(pp1);
+   Vector dshape_cx(pp1 + 1), dshape_cy(pp1 + 1);
+#endif
+   Vector dshape_ox(pp1), dshape_oy(pp1); 
+   cbasis1d.Eval(ip.x, shape_cx, dshape_cx);
+   obasis1d.Eval(ip.x, shape_ox, dshape_ox);
+   cbasis1d.Eval(ip.y, shape_cy, dshape_cy);
+   obasis1d.Eval(ip.y, shape_oy, dshape_oy);
+   gshape = 0.0; 
+   int dof2 = GetDof()/2; 
+
+   int o = 0; 
+   for (int j = 0; j < pp1; j++) 
+   {
+      for (int i = 0; i <= pp1; i++)
+      {
+         int idx, s; 
+         if ((idx = dof_map[o++]) < 0) 
+         {
+            idx = -1 - idx; 
+            s = -1; 
+         }
+         else 
+         {
+            s = +1; 
+         }
+         gshape(idx,0) = s*dshape_cx(i)*shape_oy(j); 
+         gshape(idx,1) = s*shape_cx(i)*dshape_oy(j); 
+      }
+   }
+   for (int j = 0; j <= pp1; j++) 
+   {
+		for (int i = 0; i < pp1; i++) 
+		{
+			int idx, s; 
+			if ((idx = dof_map[o++]) < 0)
+			{
+				idx = -1 - idx, s = -1; 
+			}
+			else
+			{
+				s = +1; 
+			}
+			gshape(idx,2) = s*dshape_ox(i)*shape_cy(j); 
+			gshape(idx,3) = s*shape_ox(i)*dshape_cy(j); 
+		}
+   }
+}
+
 
 const double RT_HexahedronElement::nk[18] =
 { 0.,0.,-1.,  0.,-1.,0.,  1.,0.,0.,  0.,1.,0.,  -1.,0.,0.,  0.,0.,1. };

fem/fe.hpp

@@ -2705,6 +2705,8 @@ public:
    { CalcVShape_RT(Trans, shape); }
    virtual void CalcDivShape(const IntegrationPoint &ip,
                              Vector &divshape) const;
+   virtual void CalcDShape(const IntegrationPoint &ip, 
+                             DenseMatrix &gshape) const; 
    virtual void GetLocalInterpolation(ElementTransformation &Trans,
                                       DenseMatrix &I) const
    { LocalInterpolation_RT(*this, nk, dof2nk, Trans, I); }

fem/linearform.cpp

@@ -34,6 +34,8 @@ LinearForm::LinearForm(FiniteElementSpace *f, LinearForm *lf)
 
    flfi = lf->flfi;
    flfi_marker = lf->flfi_marker;
+
+   iflfi = lf->iflfi; 
 }
 
 void LinearForm::AddDomainIntegrator(LinearFormIntegrator *lfi)
@@ -76,6 +78,11 @@ void LinearForm::AddBdrFaceIntegrator(LinearFormIntegrator *lfi,
    flfi_marker.Append(&bdr_attr_marker);
 }
 
+void LinearForm::AddInteriorFaceIntegrator(LinearFormIntegrator *lfi) 
+{
+   iflfi.Append(lfi); 
+}
+
 void LinearForm::Assemble()
 {
    Array<int> vdofs;
@@ -194,6 +201,31 @@ void LinearForm::Assemble()
          }
       }
    }
+   if (iflfi.Size()) 
+   {
+      FaceElementTransformations *tr; 
+      Mesh *mesh = fes->GetMesh(); 
+      Array<int> vdofs2; 
+      int nfaces = mesh->GetNumFaces(); 
+      for (int i = 0; i < nfaces; i++)
+      {
+         tr = mesh -> GetInteriorFaceTransformations(i); 
+         if (tr != NULL) 
+         {
+            fes -> GetElementVDofs (tr -> Elem1No, vdofs);
+            fes -> GetElementVDofs (tr -> Elem2No, vdofs2);
+            vdofs.Append (vdofs2);
+            for (int k = 0; k < iflfi.Size(); k++) 
+            {
+               iflfi[k] -> AssembleRHSElementVect(*fes -> GetFE(tr->Elem1No), 
+                                                  *fes -> GetFE(tr->Elem2No), 
+                                                  *tr, 
+                                                  elemvect); 
+               AddElementVector(vdofs, elemvect); 
+            }
+         }
+      }
+   }
 }
 
 void LinearForm::Update(FiniteElementSpace *f, Vector &v, int v_offset)
@@ -274,6 +306,7 @@ LinearForm::~LinearForm()
       for (k=0; k < dlfi.Size(); k++) { delete dlfi[k]; }
       for (k=0; k < blfi.Size(); k++) { delete blfi[k]; }
       for (k=0; k < flfi.Size(); k++) { delete flfi[k]; }
+      for (k=0; k < iflfi.Size(); k++) { delete iflfi[k]; }
    }
 }
 

fem/linearform.hpp

@@ -44,6 +44,9 @@ protected:
    Array<LinearFormIntegrator*> flfi;
    Array<Array<int>*>           flfi_marker; ///< Entries are not owned.
 
+   /// Set of interior Face Integrators to be applied. 
+   Array<LinearFormIntegrator*> iflfi; 
+
    /// The element ids where the centers of the delta functions lie
    Array<int> dlfi_delta_elem_id;
 
@@ -132,6 +135,8 @@ public:
    void AddBdrFaceIntegrator(LinearFormIntegrator *lfi,
                              Array<int> &bdr_attr_marker);
 
+   void AddInteriorFaceIntegrator(LinearFormIntegrator *lfi); 
+
    /** @brief Access all integrators added with AddDomainIntegrator() which are
        not DeltaLFIntegrator%s or they are DeltaLFIntegrator%s with non-delta
        coefficients. */
@@ -153,7 +158,7 @@ public:
    Array<Array<int>*> *GetFLFI_Marker() { return &flfi_marker; }
 
    /// Assembles the linear form i.e. sums over all domain/bdr integrators.
-   void Assemble();
+   virtual void Assemble();
 
    /// Assembles delta functions of the linear form
    void AssembleDelta();

fem/lininteg.cpp

@@ -22,6 +22,11 @@ void LinearFormIntegrator::AssembleRHSElementVect(
    mfem_error("LinearFormIntegrator::AssembleRHSElementVect(...)");
 }
 
+void LinearFormIntegrator::AssembleRHSElementVect(
+   const FiniteElement &el1, const FiniteElement &el2, FaceElementTransformations &Tr, Vector &elvect)
+{
+   mfem_error("LinearFormIntegrator::AssembleRHSElementVect(...)");
+}
 
 void DomainLFIntegrator::AssembleRHSElementVect(const FiniteElement &el,
                                                 ElementTransformation &Tr,
@@ -206,6 +211,40 @@ void BoundaryNormalLFIntegrator::AssembleRHSElementVect(
    }
 }
 
+void BoundaryNormalLFIntegrator::AssembleRHSElementVect(
+   const FiniteElement &el, FaceElementTransformations &Tr, Vector &elvect)
+{
+   int dim = el.GetDim();
+   int dof = el.GetDof();
+   Vector nor(dim), Qvec;
+
+   shape.SetSize(dof);
+   elvect.SetSize(dof);
+   elvect = 0.0;
+
+   const IntegrationRule *ir = IntRule;
+   if (ir == NULL)
+   {
+      int intorder = oa * el.GetOrder() + ob;  // <----------
+      ir = &IntRules.Get(Tr.GetGeometryType(), intorder);
+   }
+
+   for (int i = 0; i < ir->GetNPoints(); i++)
+   {
+      const IntegrationPoint &ip = ir->IntPoint(i);
+      Tr.SetAllIntPoints(&ip);
+      const IntegrationPoint &eip = Tr.GetElement1IntPoint(); 
+
+      CalcOrtho(Tr.Jacobian(), nor);
+      Q.Eval(Qvec, Tr, ip);
+
+      el.CalcShape(eip, shape);
+
+      elvect.Add(ip.weight*(Qvec*nor), shape);
+   }
+}
+
+
 void BoundaryTangentialLFIntegrator::AssembleRHSElementVect(
    const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
 {

fem/lininteg.hpp

@@ -35,6 +35,10 @@ public:
    virtual void AssembleRHSElementVect(const FiniteElement &el,
                                        FaceElementTransformations &Tr,
                                        Vector &elvect);
+   virtual void AssembleRHSElementVect(const FiniteElement &el1, 
+                                       const FiniteElement &el2, 
+                                       FaceElementTransformations &Tr, 
+                                       Vector &elvect); 
 
    virtual void SetIntRule(const IntegrationRule *ir) { IntRule = ir; }
    const IntegrationRule* GetIntRule() { return IntRule; }
@@ -182,8 +186,9 @@ public:
    virtual void AssembleRHSElementVect(const FiniteElement &el,
                                        ElementTransformation &Tr,
                                        Vector &elvect);
-
-   using LinearFormIntegrator::AssembleRHSElementVect;
+   virtual void AssembleRHSElementVect(const FiniteElement &el, 
+                                       FaceElementTransformations &Tr, 
+                                       Vector &elvec); 
 };
 
 /// Class for boundary integration \f$ L(v) = (g \cdot \tau, v) \f$ in 2D

fem/pbilinearform.cpp

@@ -442,28 +442,164 @@ void ParBilinearForm::Update(FiniteElementSpace *nfes)
    p_mat_e.Clear();
 }
 
+void ParMixedBilinearForm::pAllocMat() 
+{
+   int tr_nbr_size = trial_pfes->GetFaceNbrVSize(); 
+   int te_nbr_size = test_pfes->GetFaceNbrVSize(); 
+   if (keep_nbr_block)
+   {
+      mat = new SparseMatrix(height + te_nbr_size, width + tr_nbr_size);
+   }
+   else
+   {
+      mat = new SparseMatrix(height, width + tr_nbr_size);
+   }
+   return;
+}
+
+void ParMixedBilinearForm::AssembleSharedFaces(int skip_zeros)
+{
+   ParMesh *pmesh = trial_pfes->GetParMesh();
+   FaceElementTransformations *T;
+   Array<int> tr_vdofs1, tr_vdofs2, tr_vdofs_all;
+   Array<int> te_vdofs1, te_vdofs2, te_vdofs_all; 
+   DenseMatrix elemmat;
+
+   int nfaces = pmesh->GetNSharedFaces();
+   for (int i = 0; i < nfaces; i++)
+   {
+      T = pmesh->GetSharedFaceTransformations(i);
+      int Elem2NbrNo = T->Elem2No - pmesh->GetNE();
+      trial_pfes->GetElementVDofs(T->Elem1No, tr_vdofs1);
+      trial_pfes->GetFaceNbrElementVDofs(Elem2NbrNo, tr_vdofs2);
+      test_pfes->GetElementVDofs(T->Elem1No, te_vdofs1); 
+      test_pfes->GetFaceNbrElementVDofs(Elem2NbrNo, te_vdofs2); 
+      tr_vdofs1.Copy(tr_vdofs_all);
+      te_vdofs1.Copy(te_vdofs_all); 
+      for (int j = 0; j < tr_vdofs2.Size(); j++)
+      {
+         if (tr_vdofs2[j] >= 0)
+         {
+            tr_vdofs2[j] += width;
+         }
+         else
+         {
+            tr_vdofs2[j] -= width;
+         }
+      }
+      tr_vdofs_all.Append(tr_vdofs2);
+
+      for (int j = 0; j < te_vdofs2.Size(); j++)
+      {
+         if (te_vdofs2[j] >= 0)
+         {
+            te_vdofs2[j] += height;
+         }
+         else
+         {
+            te_vdofs2[j] -= height;
+         }
+      }
+      te_vdofs_all.Append(te_vdofs2); 
+
+      for (int k = 0; k < fbfi.Size(); k++)
+      {
+         fbfi[k]->AssembleFaceMatrix(*trial_pfes->GetFE(T->Elem1No), 
+                                     *trial_pfes->GetFE(Elem2NbrNo), 
+                                     *test_pfes->GetFE(T->Elem1No), 
+                                     *test_pfes->GetFE(Elem2NbrNo),
+                                     *T, elemmat); 
+         if (keep_nbr_block)
+         {
+            mat->AddSubMatrix(te_vdofs_all, tr_vdofs_all, elemmat, skip_zeros);
+         }
+         else
+         {
+            mat->AddSubMatrix(te_vdofs1, tr_vdofs_all, elemmat, skip_zeros);
+         }
+      }
+   }
+}
+
+void ParMixedBilinearForm::Assemble(int skip_zeros) 
+{
+   if (mat == NULL && fbfi.Size() > 0)
+   {
+      trial_pfes->ExchangeFaceNbrData();
+      test_pfes->ExchangeFaceNbrData(); 
+      pAllocMat(); 
+   }
+
+   MixedBilinearForm::Assemble(skip_zeros);
+
+   if (!ext && fbfi.Size() > 0)
+   {
+      AssembleSharedFaces(skip_zeros);
+   }
+}
 
 HypreParMatrix *ParMixedBilinearForm::ParallelAssemble()
 {
    // construct the block-diagonal matrix A
-   HypreParMatrix *A =
-      new HypreParMatrix(trial_pfes->GetComm(),
-                         test_pfes->GlobalVSize(),
-                         trial_pfes->GlobalVSize(),
-                         test_pfes->GetDofOffsets(),
-                         trial_pfes->GetDofOffsets(),
-                         mat);
+   if (fbfi.Size()==0) 
+   {
+      HypreParMatrix *A =
+         new HypreParMatrix(trial_pfes->GetComm(),
+                            test_pfes->GlobalVSize(),
+                            trial_pfes->GlobalVSize(),
+                            test_pfes->GetDofOffsets(),
+                            trial_pfes->GetDofOffsets(),
+                            mat);
 
-   HypreParMatrix *rap = RAP(test_pfes->Dof_TrueDof_Matrix(), A,
-                             trial_pfes->Dof_TrueDof_Matrix());
+      HypreParMatrix *rap = RAP(test_pfes->Dof_TrueDof_Matrix(), A,
+                                trial_pfes->Dof_TrueDof_Matrix());
 
-   delete A;
+      delete A;
 
-   return rap;
+      return rap;      
+   } 
+
+   else 
+   {
+     // handle the case when 'a' contains off-diagonal
+     int tr_lvsize = trial_pfes->GetVSize();
+     int te_lvsize = test_pfes->GetVSize(); 
+     const HYPRE_Int *tr_face_nbr_glob_ldof = trial_pfes->GetFaceNbrGlobalDofMap();
+     HYPRE_Int tr_ldof_offset = trial_pfes->GetMyDofOffset();
+
+     Array<HYPRE_Int> glob_J(mat->NumNonZeroElems());
+     int *J = mat->GetJ();
+     for (int i = 0; i < glob_J.Size(); i++)
+     {
+        if (J[i] < tr_lvsize)
+        {
+           glob_J[i] = J[i] + tr_ldof_offset;
+        }
+        else
+        {
+           glob_J[i] = tr_face_nbr_glob_ldof[J[i] - tr_lvsize];
+        }
+     }
+
+     // TODO - construct dA directly in the A format
+     HypreParMatrix *A = 
+        new HypreParMatrix(trial_pfes->GetComm(), te_lvsize, test_pfes->GlobalVSize(),
+                           trial_pfes->GlobalVSize(), mat->GetI(), glob_J,
+                           mat->GetData(), test_pfes->GetDofOffsets(),
+                           trial_pfes->GetDofOffsets());
+
+      HypreParMatrix *rap = RAP(test_pfes->Dof_TrueDof_Matrix(), A, 
+         trial_pfes->Dof_TrueDof_Matrix()); 
+
+     glob_J.DeleteAll();
+     delete A; 
+     return rap; 
+   }
 }
 
 void ParMixedBilinearForm::ParallelAssemble(OperatorHandle &A)
 {
+   if (fbfi.Size() > 0) MFEM_ABORT("parallel assemble not implemented yet for fbfi"); 
    // construct the rectangular block-diagonal matrix dA
    OperatorHandle dA(A.Type());
    dA.MakeRectangularBlockDiag(trial_pfes->GetComm(),

fem/pbilinearform.hpp

@@ -198,6 +198,10 @@ protected:
    /// Matrix and eliminated matrix
    OperatorHandle p_mat, p_mat_e;
 
+   bool keep_nbr_block = false; 
+   void pAllocMat(); 
+   void AssembleSharedFaces(int skip_zeros = 1); 
+
 private:
    /// Copy construction is not supported; body is undefined.
    ParMixedBilinearForm(const ParMixedBilinearForm &);
@@ -238,6 +242,11 @@ public:
       test_pfes  = test_fes;
    }
 
+   void KeepNbrBlock(bool knb = true) { keep_nbr_block = knb; }
+
+   /// Assemble the local matrix
+   void Assemble(int skip_zeros = 1); 
+
    /// Returns the matrix assembled on the true dofs, i.e. P_test^t A P_trial.
    HypreParMatrix *ParallelAssemble();
 

fem/pgridfunc.cpp

@@ -471,6 +471,46 @@ void ParGridFunction::GetVectorValue(ElementTransformation &T,
    }
 }
 
+void ParGridFunction::GetGradient(ElementTransformation &T, Vector &grad) const 
+{
+   Array<int> dofs;
+   Vector DofVal, LocVec;
+   int nbr_el_no = T.ElementNo - pfes->GetParMesh()->GetNE();
+   if (nbr_el_no >= 0) 
+   {
+      switch (T.ElementType)
+      {
+         case ElementTransformation::ELEMENT:
+         {
+            pfes->GetFaceNbrElementVDofs(nbr_el_no, dofs);
+            const FiniteElement *fe = pfes->GetFaceNbrFE(nbr_el_no);
+            MFEM_ASSERT(fe->GetMapType() == FiniteElement::VALUE,
+                        "invalid FE map type");
+            int spaceDim = pfes->GetMesh()->SpaceDimension();
+            int dim = fe->GetDim(), dof = fe->GetDof();
+            DenseMatrix dshape(dof, dim);
+            Vector lval, gh(dim);
+
+            grad.SetSize(spaceDim);
+            face_nbr_data.GetSubVector(dofs, LocVec);
+            fe->CalcDShape(T.GetIntPoint(), dshape);
+            dshape.MultTranspose(LocVec, gh);
+            T.InverseJacobian().MultTranspose(gh, grad);
+         }
+         break;
+         default:
+         {
+            MFEM_ABORT("GridFunction::GetGradient: Unsupported element type \""
+                       << T.ElementType << "\"");
+         }
+      }
+   }
+   else 
+   {
+      GridFunction::GetGradient(T, grad); 
+   }
+}
+
 void ParGridFunction::ProjectCoefficient(Coefficient &coeff)
 {
    DeltaCoefficient *delta_c = dynamic_cast<DeltaCoefficient *>(&coeff);

fem/pgridfunc.hpp

@@ -221,6 +221,8 @@ public:
                                const IntegrationPoint &ip,
                                Vector &val, Vector *tr = NULL) const;
 
+   virtual void GetGradient(ElementTransformation &T, Vector &grad) const; 
+
    using GridFunction::ProjectCoefficient;
    virtual void ProjectCoefficient(Coefficient &coeff);
 

fem/plinearform.cpp

@@ -43,6 +43,58 @@ void ParLinearForm::MakeRef(ParFiniteElementSpace *pf, Vector &v, int v_offset)
    pfes = pf;
 }
 
+void ParLinearForm::AssembleSharedFaces() 
+{
+	ParMesh *pmesh = pfes->GetParMesh();
+	FaceElementTransformations *T;
+	Array<int> vdofs1, vdofs2, vdofs_all;
+	Vector elvec;
+
+	int nfaces = pmesh->GetNSharedFaces();
+	for (int i = 0; i < nfaces; i++)
+	{
+	   T = pmesh->GetSharedFaceTransformations(i);
+	   int Elem2NbrNo = T->Elem2No - pmesh->GetNE();
+	   pfes->GetElementVDofs(T->Elem1No, vdofs1);
+	   pfes->GetFaceNbrElementVDofs(Elem2NbrNo, vdofs2);
+	   vdofs1.Copy(vdofs_all);
+	   int height = pfes->GetVSize(); 
+	   for (int j = 0; j < vdofs2.Size(); j++)
+	   {
+	      if (vdofs2[j] >= 0)
+	      {
+	         vdofs2[j] += height;
+	      }
+	      else
+	      {
+	         vdofs2[j] -= height;
+	      }
+	   }
+	   vdofs_all.Append(vdofs2);
+	   for (int k = 0; k < iflfi.Size(); k++)
+	   {
+	      iflfi[k]->AssembleRHSElementVect(*pfes->GetFE(T->Elem1No),
+	      	*pfes->GetFaceNbrFE(Elem2NbrNo),
+	      	*T, elvec); 
+	      Vector local; 
+	      local.MakeRef(elvec, 0, vdofs1.Size()); 
+	      AddElementVector(vdofs1, local); 
+	   }
+	}
+
+}
+
+void ParLinearForm::Assemble() 
+{
+	if (iflfi.Size()>0) {
+		pfes->ExchangeFaceNbrData(); 
+	}
+	LinearForm::Assemble(); 
+	if (iflfi.Size()>0) {
+		AssembleSharedFaces(); 
+	}
+}
+
 void ParLinearForm::ParallelAssemble(Vector &tv)
 {
    const Operator* prolong = pfes->GetProlongationMatrix();

fem/plinearform.hpp

@@ -27,6 +27,7 @@ class ParLinearForm : public LinearForm
 {
 protected:
    ParFiniteElementSpace *pfes; ///< Points to the same object as #fes
+   void AssembleSharedFaces(); 
 
 private:
    /// Copy construction is not supported; body is undefined.
@@ -113,6 +114,8 @@ public:
        build option MFEM_DEBUG is enabled. */
    void MakeRef(ParFiniteElementSpace *pf, Vector &v, int v_offset);
 
+   void Assemble(); 
+
    /// Assemble the vector on the true dofs, i.e. P^t v.
    void ParallelAssemble(Vector &tv);
 

linalg/blockoperator.cpp

@@ -381,4 +381,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

@@ -226,7 +226,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).
    /**
@@ -282,6 +282,180 @@ 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

@@ -128,6 +128,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);
@@ -915,17 +922,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
@@ -966,6 +982,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)
 {
@@ -1568,6 +1624,40 @@ 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;
+}
+
 #if MFEM_HYPRE_VERSION < 21400
 
 HypreParMatrix *Add(double alpha, const HypreParMatrix &A,
@@ -1608,6 +1698,16 @@ HypreParMatrix *Add(double alpha, const HypreParMatrix &A,
    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);
+
+   return new HypreParMatrix(C_hypre);
+}
+
 HypreParMatrix * ParAdd(const HypreParMatrix *A, const HypreParMatrix *B)
 {
    hypre_ParCSRMatrix *C;
@@ -2184,9 +2284,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;
 }
 
@@ -2490,6 +2590,8 @@ HypreSolver::HypreSolver()
 {
    A = NULL;
    setup_called = 0;
+   final_res_norm = -1;
+   num_iterations = -1;
    B = X = NULL;
    error_mode = ABORT_HYPRE_ERRORS;
 }
@@ -2499,6 +2601,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;
 }
@@ -2669,8 +2773,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;
 
@@ -2715,6 +2817,9 @@ void HyprePCG::Mult(const HypreParVector &b, HypreParVector &x) const
    b.HostRead();
    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)
    {
@@ -2726,10 +2831,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)
@@ -2804,6 +2905,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);
@@ -2838,8 +2945,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;
 
@@ -2879,6 +2984,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)
    {
@@ -2887,10 +2995,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)
@@ -3339,6 +3443,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:
@@ -3607,6 +3777,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++)
@@ -4674,4 +5088,4 @@ HypreAME::StealEigenvectors()
 
 }
 
-#endif
+#endif

linalg/hypre.hpp

@@ -81,6 +81,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
@@ -113,9 +116,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
@@ -228,22 +231,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
@@ -328,6 +334,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; }
@@ -387,6 +394,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. */
@@ -397,6 +406,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. */
@@ -650,7 +670,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();
@@ -661,7 +681,7 @@ public:
                  double poly_fraction = .3, int eig_est_cg_iter = 10);
 
    /// 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
@@ -717,6 +737,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;
@@ -726,6 +748,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;
 
@@ -755,6 +780,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
 {
@@ -832,6 +876,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);
@@ -1116,9 +1161,83 @@ public:
        construct A. */
    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; }
 
@@ -1127,6 +1246,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

@@ -1225,7 +1225,8 @@ void BiCGSTABSolver::Mult(const Vector &b, Vector &x) const
       final_norm = resid;
       final_iter = 0;
       converged = 1;
-      return;
+      // return;
+      goto finish; 
    }
 
    for (i = 1; i <= max_iter; i++)
@@ -1242,7 +1243,8 @@ void BiCGSTABSolver::Mult(const Vector &b, Vector &x) const
          final_norm = resid;
          final_iter = i;
          converged = 0;
-         return;
+         // return;
+         goto finish; 
       }
       if (i == 1)
       {
@@ -1276,7 +1278,8 @@ void BiCGSTABSolver::Mult(const Vector &b, Vector &x) const
          final_norm = resid;
          final_iter = i;
          converged = 1;
-         return;
+         // return;
+         goto finish; 
       }
       if (print_level >= 0)
          mfem::out << "   Iteration : " << setw(3) << i
@@ -1309,20 +1312,38 @@ void BiCGSTABSolver::Mult(const Vector &b, Vector &x) const
          final_norm = resid;
          final_iter = i;
          converged = 1;
-         return;
+         // return;
+         goto finish; 
       }
       if (omega == 0)
       {
          final_norm = resid;
          final_iter = i;
          converged = 0;
-         return;
+         // return;
+         goto finish; 
       }
    }
 
    final_norm = resid;
    final_iter = max_iter;
    converged = 0;
+finish:
+   if (print_level == 1 || print_level == 3)
+   {
+      mfem::out << "   Iteration : " << setw(3) << final_iter
+                << "  ||B r|| = " << final_norm << '\n';
+   }
+   else if (print_level == 2)
+   {
+      mfem::out << "BiCG: Number of iterations: " << final_iter << '\n';
+   }
+   if (print_level >= 0 && !converged)
+   {
+      mfem::out << "BiCG: No convergence!\n";
+   }
+
+   Monitor(final_iter, final_norm, r, x, true);
 }
 
 int BiCGSTAB(const Operator &A, Vector &x, const Vector &b, Solver &M,

linalg/vector.cpp

@@ -522,6 +522,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();
@@ -588,6 +604,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

@@ -302,6 +302,10 @@ public:
        the -val in from this Vector. */
    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 @a dofs to the given @a value.
    /** Negative dof values cause the -dof-1 position in this Vector to receive
        the -value. */
@@ -317,6 +321,10 @@ public:
        Vector to receive the -val from @a elem_data. */
    void SetSubVector(const Array<int> &dofs, double *elem_data);
 
+   // ADDED //
+   void SetSubVector(int index_low, int index_high, Vector &elemvect);
+   // ADDED //
+
    /** @brief Add elements of the @a elemvect Vector to the entries listed in @a
        dofs. Negative dof values cause the -dof-1 position in this Vector to add
        the -val from @a elemvect. */
@@ -333,6 +341,10 @@ public:
    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 @a dofs Array to the given @a val.
    void SetSubVectorComplement(const Array<int> &dofs, const double val);
 

mesh/mesh.cpp

@@ -965,6 +965,94 @@ FaceElementTransformations *Mesh::GetFaceElementTransformations(int FaceNo,
    return &FaceElemTr;
 }
 
+FaceElementTransformations* Mesh::GetFaceElementTransformations(FaceElementTransformationsData &fetd, int FaceNo, int mask) {
+   FaceInfo &face_info = faces_info[FaceNo];
+   FaceElementTransformations &face = fetd.face; 
+
+   int cmask = 0;
+   face.SetConfigurationMask(cmask);
+   face.Elem1 = NULL;
+   face.Elem2 = NULL;
+
+   // setup the transformation for the first element
+   face.Elem1No = face_info.Elem1No;
+   if (mask & FaceElementTransformations::HAVE_ELEM1)
+   {
+      GetElementTransformation(face.Elem1No, &fetd.Elem1);
+      face.Elem1 = &fetd.Elem1;
+      cmask |= 1;
+   }
+
+   //  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"
+   face.Elem2No = face_info.Elem2No;
+   if ((mask & FaceElementTransformations::HAVE_ELEM2) &&
+       face.Elem2No >= 0)
+   {
+#ifdef MFEM_DEBUG
+      if (NURBSext && (mask & FaceElementTransformations::HAVE_ELEM1))
+      { MFEM_ABORT("NURBS mesh not supported!"); }
+#endif
+      GetElementTransformation(face.Elem2No, &fetd.Elem2);
+      face.Elem2 = &fetd.Elem2;
+      cmask |= 2;
+   }
+
+   // setup the face transformation
+   if (mask & FaceElementTransformations::HAVE_FACE)
+   {
+      GetFaceTransformation(FaceNo, &face);
+      cmask |= 16;
+   }
+   else
+   {
+      face.SetGeometryType(GetFaceGeometryType(FaceNo));
+   }
+
+   // setup Loc1 & Loc2
+   int face_type = GetFaceElementType(FaceNo);
+   if (mask & FaceElementTransformations::HAVE_LOC1)
+   {
+      int elem_type = GetElementType(face_info.Elem1No);
+      GetLocalFaceTransformation(face_type, elem_type,
+                                 face.Loc1.Transf, face_info.Elem1Inf);
+      cmask |= 4;
+   }
+   if ((mask & FaceElementTransformations::HAVE_LOC2) &&
+       face.Elem2No >= 0)
+   {
+      int elem_type = GetElementType(face_info.Elem2No);
+      GetLocalFaceTransformation(face_type, elem_type,
+                                 face.Loc2.Transf, face_info.Elem2Inf);
+
+      // NC meshes: prepend slave edge/face transformation to Loc2
+      if (Nonconforming() && IsSlaveFace(face_info))
+      {
+         ApplyLocalSlaveTransformation(face, face_info, false);
+      }
+      cmask |= 8;
+   }
+
+   face.SetConfigurationMask(cmask);
+
+   // This check can be useful for internal debugging, however it will fail on
+   // periodic boundary faces, so we keep it disabled in general.
+#if 0
+#ifdef MFEM_DEBUG
+   double dist = FaceElemTr.CheckConsistency();
+   if (dist >= 1e-12)
+   {
+      mfem::out << "\nInternal error: face id = " << FaceNo
+                << ", dist = " << dist << '\n';
+      FaceElemTr.CheckConsistency(1); // print coordinates
+      MFEM_ABORT("internal error");
+   }
+#endif
+#endif
+   return &fetd; 
+}
+
 bool Mesh::IsSlaveFace(const FaceInfo &fi) const
 {
    return fi.NCFace >= 0 && nc_faces_info[fi.NCFace].Slave;
@@ -1037,6 +1125,39 @@ FaceElementTransformations *Mesh::GetBdrFaceTransformations(int BdrElemNo)
    return tr;
 }
 
+FaceElementTransformations* Mesh::GetInteriorFaceTransformations(FaceElementTransformationsData &fetd, int FaceNo) 
+{
+	if (faces_info[FaceNo].Elem2No < 0) return NULL; 
+	else return GetFaceElementTransformations(fetd, FaceNo); 
+}
+
+FaceElementTransformations* Mesh::GetBdrFaceTransformations(FaceElementTransformationsData &fetd, int BdrElemNo)
+{
+	int fn;
+	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)
+	{
+		return NULL; 
+	}
+	GetFaceElementTransformations(fetd, fn, 21);
+	fetd.face.Attribute = boundary[BdrElemNo]->GetAttribute();
+	fetd.face.ElementNo = BdrElemNo;
+	fetd.face.ElementType = ElementTransformation::BDR_FACE;
+	return &fetd;  
+}
+
 void Mesh::GetFaceElements(int Face, int *Elem1, int *Elem2) const
 {
    *Elem1 = faces_info[Face].Elem1No;

mesh/mesh.hpp

@@ -40,6 +40,7 @@ class NURBSExtension;
 class FiniteElementSpace;
 class GridFunction;
 struct Refinement;
+class FaceElementTransformationsData; 
 
 /** An enum type to specify if interior or boundary faces are desired. */
 enum class FaceType : bool {Interior, Boundary};
@@ -1029,6 +1030,12 @@ public:
 
    FaceElementTransformations *GetBdrFaceTransformations (int BdrElemNo);
 
+   // --- added --- 
+   FaceElementTransformations* GetFaceElementTransformations(FaceElementTransformationsData &fetd, int FaceNo, int mask=31); 
+   FaceElementTransformations* GetInteriorFaceTransformations(FaceElementTransformationsData &fetd, int FaceNo); 
+   FaceElementTransformations* GetBdrFaceTransformations(FaceElementTransformationsData &fetd, int BdrElemNo); 
+   // --- added --- 
+
    /// Return true if the given face is interior. @sa FaceIsTrueInterior().
    bool FaceIsInterior(int FaceNo) const
    {
@@ -1520,6 +1527,13 @@ inline void ShiftRight(int &a, int &b, int &c)
    a = c;  c = b;  b = t;
 }
 
+class FaceElementTransformationsData {
+public:
+	IsoparametricTransformation Elem1, Elem2; 
+	FaceElementTransformations face; 
+	FaceElementTransformations* operator&() { return &face; }
+}; 
+
 }
 
 #endif

mesh/pmesh.cpp

@@ -318,12 +318,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 //
       }
    }
 
@@ -2510,6 +2518,116 @@ GetSharedFaceTransformations(int sf, bool fill2)
    return &FaceElemTr;
 }
 
+FaceElementTransformations* ParMesh::GetSharedFaceTransformations(FaceElementTransformationsData &fetd, int sf, bool fill2)
+{
+   int FaceNo = GetSharedFace(sf);
+
+   FaceInfo &face_info = faces_info[FaceNo];
+
+   bool is_slave = Nonconforming() && IsSlaveFace(face_info);
+   bool is_ghost = Nonconforming() && FaceNo >= GetNumFaces();
+
+   FaceElementTransformations &face = fetd.face; 
+
+   int mask = 0;
+   face.SetConfigurationMask(0);
+   face.Elem1 = NULL;
+   face.Elem2 = NULL;
+
+   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
+   face.Elem1No = face_info.Elem1No;
+   GetElementTransformation(face.Elem1No, &fetd.Elem1);
+   face.Elem1 = &fetd.Elem1;
+   mask |= FaceElementTransformations::HAVE_ELEM1;
+
+   // setup the transformation for the second (neighbor) element
+   int Elem2NbrNo;
+   if (fill2)
+   {
+      Elem2NbrNo = -1 - face_info.Elem2No;
+      // Store the "shifted index" for element 2 in FaceElemTr.Elem2No.
+      // `Elem2NbrNo` is the index of the face neighbor (starting from 0),
+      // and `FaceElemTr.Elem2No` will be offset by the number of (local)
+      // elements in the mesh.
+      face.Elem2No = NumOfElements + Elem2NbrNo;
+      GetFaceNbrElementTransformation(Elem2NbrNo, &fetd.Elem2);
+      face.Elem2 = &fetd.Elem2;
+      mask |= FaceElementTransformations::HAVE_ELEM2;
+   }
+   else
+   {
+      face.Elem2No = -1;
+   }
+
+   // setup the face transformation if the face is not a ghost
+   if (!is_ghost)
+   {
+      GetFaceTransformation(FaceNo, &face);
+      // NOTE: The above call overwrites FaceElemTr.Loc1
+      mask |= FaceElementTransformations::HAVE_FACE;
+   }
+   else
+   {
+      face.SetGeometryType(face_geom);
+   }
+
+   // setup Loc1 & Loc2
+   int elem_type = GetElementType(face_info.Elem1No);
+   GetLocalFaceTransformation(face_type, elem_type, face.Loc1.Transf,
+                              face_info.Elem1Inf);
+   mask |= FaceElementTransformations::HAVE_LOC1;
+
+   if (fill2)
+   {
+      elem_type = face_nbr_elements[Elem2NbrNo]->GetType();
+      GetLocalFaceTransformation(face_type, elem_type, face.Loc2.Transf,
+                                 face_info.Elem2Inf);
+      mask |= FaceElementTransformations::HAVE_LOC2;
+   }
+
+   // adjust Loc1 or Loc2 of the master face if this is a slave face
+   if (is_slave)
+   {
+      if (is_ghost || fill2)
+      {
+         // is_ghost -> modify side 1, otherwise -> modify side 2:
+         ApplyLocalSlaveTransformation(face, face_info, is_ghost);
+      }
+   }
+
+   // for ghost faces we need a special version of GetFaceTransformation
+   if (is_ghost)
+   {
+      GetGhostFaceTransformation(&face, face_type, face_geom);
+      mask |= FaceElementTransformations::HAVE_FACE;
+   }
+
+   face.SetConfigurationMask(mask);
+
+   // This check can be useful for internal debugging, however it will fail on
+   // periodic boundary faces, so we keep it disabled in general.
+#if 0
+#ifdef MFEM_DEBUG
+   double dist = FaceElemTr.CheckConsistency();
+   if (dist >= 1e-12)
+   {
+      mfem::out << "\nInternal error: face id = " << FaceNo
+                << ", dist = " << dist << ", rank = " << MyRank << '\n';
+      FaceElemTr.CheckConsistency(1); // print coordinates
+      MFEM_ABORT("internal error");
+   }
+#endif
+#endif
+   return &fetd; 
+}
+
 int ParMesh::GetNSharedFaces() const
 {
    if (Conforming())

mesh/pmesh.hpp

@@ -257,10 +257,43 @@ 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; 
+       }
+   }
+
    ///@{ @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); }
@@ -298,6 +331,10 @@ public:
    FaceElementTransformations *
    GetSharedFaceTransformations(int sf, bool fill2 = true);
 
+   // --- added --- 
+   FaceElementTransformations* GetSharedFaceTransformations(FaceElementTransformationsData &fetd, int sf, bool fill2 = true); 
+   // --- added --- 
+
    ElementTransformation *
    GetFaceNbrElementTransformation(int i)
    {