WIP: updates to MG setup

for EA capabilities

CMakeLists.txt

@@ -15,7 +15,7 @@ set(USER_CONFIG "${CMAKE_CURRENT_SOURCE_DIR}/config/user.cmake" CACHE PATH
   "Path to optional user configuration file.")
 
 # Require C++11 and disable compiler-specific extensions
-set(CMAKE_CXX_STANDARD 11)
+set(CMAKE_CXX_STANDARD 14)
 set(CMAKE_CXX_STANDARD_REQUIRED ON)
 set(CMAKE_CXX_EXTENSIONS OFF)
 
@@ -95,7 +95,7 @@ if (MFEM_USE_CUDA)
       set(CMAKE_CUDA_HOST_COMPILER ${CMAKE_CXX_COMPILER})
    endif()
    enable_language(CUDA)
-   set(CMAKE_CUDA_STANDARD 11)
+   set(CMAKE_CUDA_STANDARD 14)
    set(CMAKE_CUDA_STANDARD_REQUIRED ON)
    set(CMAKE_CUDA_EXTENSIONS OFF)
    set(CUDA_FLAGS "--expt-extended-lambda")

examples/ginkgo/CMakeLists.txt

@@ -12,6 +12,9 @@
 set(GINKGO_EXAMPLES_SRCS)
 list(APPEND GINKGO_EXAMPLES_SRCS
   ex1.cpp
+  ex1lor.cpp
+  ex1lorcoeff.cpp
+  ex1tetlorcoeff.cpp
 )
 
 # Include the source directory where mfem.hpp and mfem-performance.hpp are.
@@ -33,6 +36,8 @@ add_mfem_examples(GINKGO_EXAMPLES_SRCS ${PFX} "" test_ginkgo)
 
 # Command line options for the tests.
 set(EX1_COMMON_OPTS ex1 -m ../data/star.mesh --use_gko_solver)
+set(EX1LOR_COMMON_OPTS ex1lor -m ../../data/beam-hex.mesh -pc)
+
 set(EX1_TEST_OPTS ${EX9_COMMON_OPTS})
 
 # Add the tests: one test per source file.

examples/ginkgo/amgx.json

@@ -0,0 +1,17 @@
+{
+    "config_version": 2,
+    "solver": {
+        "solver": "AMG",
+        "smoother": {
+            "scope": "jacobi",
+            "solver": "JACOBI_L1"
+        },
+        "presweeps": 1,
+        "interpolator": "D2",
+        "max_row_sum" : 0.9,
+        "strength_threshold" : 0.25,
+        "postsweeps": 1,
+        "max_iters": 1,
+        "cycle": "V"
+    }
+}

examples/ginkgo/ex1lor.cpp

@@ -0,0 +1,721 @@
+//                          MFEM Example 1, modified
+//
+// This code has been modified from `ex1.cpp` provided in the examples of
+// MFEM.  The sections not marked as related to Ginkgo are largely unchanged
+// from the version provided by MFEM.
+//
+// This code also contains portions from `miniapps/performance/ex1.cpp`
+// pertaining to the LOR preconditioner.  The preconditioner types
+// used in this example are applied to the SparseMatrix created from
+// the LOR mesh, so they are all LOR preconditioners with different
+// subtypes.
+//
+// The default mesh option is "beam-hex.mesh", provided by MFEM.
+// Important non-default options:
+//  -m [file] : Mesh file.
+//  -d "cuda" : Use the MFEM cuda backend and Ginkgo CudaExecutor.
+//  -pc-type "gko:ilu" : Use the Ginkgo ILU preconditioner (default is Block
+//                       Jacobi)
+//  -pc-type "none" : No LOR preconditioner
+//
+//  Options only for the Block Jacobi preconditioner (default:)
+//  -pc-so "none" : Don't let Ginkgo automatically pick options for precision
+//                  reduction in the storage of the Block Jacobi preconditioner
+//  -pc-acc [value] : Accuracy parameter.
+//
+// MFEM's provided information about `ex1.cpp`:
+// Description:  This example code demonstrates the use of MFEM to define a
+//               simple finite element discretization of the Laplace problem
+//               -Delta u = 1 with homogeneous Dirichlet boundary conditions.
+//               Specifically, we discretize using a FE space of the specified
+//               order, or if order < 1 using an isoparametric/isogeometric
+//               space (i.e. quadratic for quadratic curvilinear mesh, NURBS for
+//               NURBS mesh, etc.)
+//
+//               The example highlights the use of mesh refinement, finite
+//               element grid functions, as well as linear and bilinear forms
+//               corresponding to the left-hand side and right-hand side of the
+//               discrete linear system. We also cover the explicit elimination
+//               of essential boundary conditions, static condensation, and the
+//               optional connection to the GLVis tool for visualization.
+
+#include "mfem.hpp"
+
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+void PermuteSparseMatrix(SparseMatrix &A, Array<int> &pinv, double shift=0.0)
+{
+   int n = pinv.Size();
+
+   SparseMatrix PA(A);
+
+   Array<int> p(n);
+   for (int i=0; i<n; ++i)
+   {
+      p[pinv[i]] = i;
+   }
+
+   // Set LU = A(P,P) using the permutation generated above
+   const int *IA = A.GetI();
+   const int *JA = A.GetJ();
+   const double *VA = A.GetData();
+
+   int *I = PA.GetI();
+   int *J = PA.GetJ();
+   double *V = PA.GetData();
+
+   I[0] = 0;
+   for (int i=0; i<n; ++i)
+   {
+      int pi = p[i];
+      int nnz_pi = IA[pi+1] - IA[pi];
+      I[i+1] = I[i] + nnz_pi;
+      for (int jj=0; jj<nnz_pi; ++jj)
+      {
+         int pj = JA[IA[pi] + jj];
+         int j = pinv[pj];
+
+         J[I[i] + jj] = j;
+         V[I[i] + jj] = VA[IA[pi] + jj];
+         if (i == j) { V[I[i] + jj] += shift; }
+      }
+   }
+
+   PA.SortColumnIndices();
+   PA.Swap(A);
+}
+
+// helper functions for cg/pcg solve with timer and iter count return
+int cg_solve(const Operator &A, const Vector &b, Vector &x,
+             int print_iter, int max_num_iter,
+             double RTOLERANCE, double ATOLERANCE, double &it_time)
+{
+
+   CGSolver cg;
+   cg.SetPrintLevel(print_iter);
+   cg.SetMaxIter(max_num_iter);
+   cg.SetRelTol(sqrt(RTOLERANCE));
+   cg.SetAbsTol(sqrt(ATOLERANCE));
+   cg.SetOperator(A);
+
+   tic_toc.Clear();
+   tic_toc.Start();
+
+   cg.Mult(b, x);
+
+   tic_toc.Stop();
+   it_time = tic_toc.RealTime();
+
+   return cg.GetNumIterations();
+}
+
+int pcg_solve(const Operator &A, Solver &B, const Vector &b, Vector &x,
+              int print_iter, int max_num_iter,
+              double RTOLERANCE, double ATOLERANCE, double &it_time)
+{
+
+   CGSolver pcg;
+   pcg.SetPrintLevel(print_iter);
+   pcg.SetMaxIter(max_num_iter);
+   pcg.SetRelTol(sqrt(RTOLERANCE));
+   pcg.SetAbsTol(sqrt(ATOLERANCE));
+   pcg.SetOperator(A);
+   pcg.SetPreconditioner(B);
+
+   tic_toc.Clear();
+   tic_toc.Start();
+
+   pcg.Mult(b, x);
+
+   tic_toc.Stop();
+   it_time = tic_toc.RealTime();
+
+   return pcg.GetNumIterations();
+}
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   const char *mesh_file = "../../data/beam-hex.mesh";
+   int ref_levels = 3;
+   int order = 2;
+   const char *basis_type = "G"; // Gauss-Lobatto
+   bool static_cond = false;
+   bool pa = true;
+   const char *device_config = "cpu";
+   bool visualization = true;
+   const char *pc_type = "gko:bj";
+   const char *pc_storage_opt = "auto";
+   double pc_acc = 1.e-1;
+   int pc_max_bs = 32;
+   int permute = 0;
+   bool output_sol = false;
+   bool output_pc = false;
+   int isai_sparsity_power = 1;
+   int par_ilu_its = 0;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
+   args.AddOption(&ref_levels, "-l", "--refinement-levels",
+                  "Number of uniform refinement levels for mesh.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
+                  "--no-static-condensation", "Enable static condensation.");
+   args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
+                  "--no-partial-assembly", "Enable Partial Assembly.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&pc_type, "-pc-type", "--preconditioner-type",
+                  "Type of preconditioner used on LOR matrix.");
+   args.AddOption(&pc_storage_opt, "-pc-so",
+                  "--preconditioner-storage-optimization",
+                  "Type of precision storage optimization to use for Ginkgo BlockJacobi.");
+   args.AddOption(&pc_acc, "-pc-acc", "--preconditioner-accuracy",
+                  "Accuracy parameter for Ginkgo BlockJacobi.");
+   args.AddOption(&pc_max_bs, "-pc-mbs", "--preconditioner-max-block-size",
+                  "Maximum block size for Ginkgo BlockJacobi.");
+   args.AddOption(&permute, "-per", "--permutation",
+                  "Specify preconditioner permutation.");
+   args.AddOption(&output_sol, "-out", "--output-solution-and-mesh", "-no-out",
+                  "--no-solution-and-mesh-output",
+                  "Output mesh and solution for inspection.");
+   args.AddOption(&output_pc, "-out-pc", "--output-lor-matrix-and-mesh",
+                  "-no-out-pc",
+                  "--no-lor-matrix-and-mesh-output",
+                  "Output LOR mesh and sparse matrix for inspection.");
+   args.AddOption(&isai_sparsity_power, "-isai-sp", "--isai-sparsity-power",
+                  "Power to use for sparsity pattern of ISAI in Ginkgo ILU-ISAI.");
+   args.AddOption(&par_ilu_its, "-pilu-its", "--par-ilu-iterations",
+                  "Number of iterations for the Ginkgo ParILU algorithm.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   enum PCType { NONE, GKO_BLOCK_JACOBI, GKO_ILU, GKO_ILU_ISAI, GKO_CUILU, GKO_CUILU_ISAI, MFEM_GS, MFEM_UMFPACK };
+   PCType pc_choice;
+   bool pc = true;
+   const char *trisolve_type = "exact"; //only used for ILU
+   if (!strcmp(pc_type, "gko:bj")) { pc_choice = GKO_BLOCK_JACOBI; }
+   else if (!strcmp(pc_type, "gko:ilu")) { pc_choice = GKO_ILU; }
+   else if (!strcmp(pc_type, "gko:ilu-isai"))
+   {
+      pc_choice = GKO_ILU_ISAI;
+      trisolve_type = "isai";
+   }
+   else if (!strcmp(pc_type, "gko:cuilu")) { pc_choice = GKO_CUILU; }
+   else if (!strcmp(pc_type, "gko:cuilu-isai"))
+   {
+      pc_choice = GKO_CUILU_ISAI;
+      trisolve_type = "isai";
+   }
+   else if (!strcmp(pc_type, "mfem:gs")) { pc_choice = MFEM_GS; }
+   else if (!strcmp(pc_type, "mfem:umf"))
+   {
+#ifdef MFEM_USE_SUITESPARSE
+      pc_choice = MFEM_UMFPACK;
+#else
+      mfem_error("Preconditioner requires SuiteSparse");
+#endif
+   }
+   else if (!strcmp(pc_type, "none"))
+   {
+      pc_choice = NONE;
+      pc = false;
+   }
+   else
+   {
+      mfem_error("Invalid Preconditioner specified");
+      return 3;
+   }
+
+   // ---------------------------------------------------------------
+   // -------------------- Start Ginkgo set-up ----------------------
+
+   // Create Ginkgo executor.
+
+   // This will point to the selected executor default executor
+   std::shared_ptr<gko::Executor> executor;
+
+   // We will always need an OpenMP executor.
+   auto omp_executor = gko::OmpExecutor::create();
+
+   // If the user has requested to use CUDA, then build a
+   // CudaExecutor and set `executor` to it; otherwise,
+   // use the OmpExecutor
+   if (!strcmp(device_config, "cuda"))
+   {
+      auto cuda_executor =
+         gko::CudaExecutor::create(0, gko::OmpExecutor::create());
+      executor = cuda_executor;
+   }
+   else
+   {
+      executor = omp_executor;
+   }
+
+   // --------------------- End Ginkgo set-up -----------------------
+   // ---------------------------------------------------------------
+
+   // See class BasisType in fem/fe_coll.hpp for available basis types
+   int basis = BasisType::GetType(basis_type[0]);
+   cout << "Using " << BasisType::Name(basis) << " basis ..." << endl;
+
+   // 2. Enable hardware devices such as GPUs, and programming models such as
+   //    CUDA, OCCA, RAJA and OpenMP based on command line options.
+   Device device(device_config);
+   device.Print();
+
+   // 3. Read the mesh from the given mesh file. We can handle triangular,
+   //    quadrilateral, tetrahedral, hexahedral, surface and volume meshes with
+   //    the same code.
+   Mesh *mesh = new Mesh(mesh_file, 1, 1);
+   int dim = mesh->Dimension();
+
+   // 4. Refine the mesh to increase the resolution. In this example we do
+   //    'ref_levels' of uniform refinement.
+   {
+      for (int l = 0; l < ref_levels; l++)
+      {
+         mesh->UniformRefinement();
+      }
+   }
+   cout << "Total elements in refined mesh: " <<  mesh->GetNE() << std::endl;
+
+   // 5. Define a finite element space on the mesh. Here we use continuous
+   //    Lagrange finite elements of the specified order. If order < 1, we
+   //    instead use an isoparametric/isogeometric space.
+   FiniteElementCollection *fec;
+   if (order > 0)
+   {
+      fec = new H1_FECollection(order, dim, basis);
+   }
+   else if (mesh->GetNodes())
+   {
+      fec = mesh->GetNodes()->OwnFEC();
+      cout << "Using isoparametric FEs: " << fec->Name() << endl;
+   }
+   else
+   {
+      fec = new H1_FECollection(order = 1, dim, basis);
+   }
+   FiniteElementSpace *fespace = new FiniteElementSpace(mesh, fec);
+   cout << "Number of finite element unknowns: " << fespace->GetTrueVSize()
+        << endl;
+
+   // Create the LOR mesh and finite element space. In the settings of this
+   // example, we can transfer between HO and LOR with the identity operator.
+   Mesh *mesh_lor = NULL;
+   FiniteElementCollection *fec_lor = NULL;
+   FiniteElementSpace *fespace_lor = NULL;
+   Array<int> *inv_reordering = NULL;
+   if (pc)
+   {
+      int basis_lor = basis;
+      if (basis == BasisType::Positive) { basis_lor=BasisType::ClosedUniform; }
+      mesh_lor = new Mesh(mesh, order, basis_lor);
+      fec_lor = new H1_FECollection(1, dim);
+      fespace_lor = new FiniteElementSpace(mesh_lor, fec_lor);
+
+      if (permute == 1)
+      {
+
+         tic_toc.Clear();
+         tic_toc.Start();
+
+         const Table &pre_reorder_dofs = fespace_lor->GetElementToDofTable();
+         const Table pre_reorder_dofs_copy(pre_reorder_dofs);
+         fespace_lor->ReorderElementToDofTable();
+         const Table &post_reorder_dofs = fespace_lor->GetElementToDofTable();
+
+         inv_reordering = new Array<int>(fespace_lor->GetTrueVSize());
+         for (int i = 0; i < pre_reorder_dofs.Size(); i++)
+         {
+
+            Array<int> old_row;
+            Array<int> new_row;
+            pre_reorder_dofs_copy.GetRow(i, old_row);
+            post_reorder_dofs.GetRow(i, new_row);
+            for (int j = 0; j < pre_reorder_dofs_copy.RowSize(i); j++)
+            {
+               int new_dof = new_row[j];
+               int old_dof = old_row[j];
+               (*inv_reordering)[old_dof] = new_dof;
+            }
+         }
+
+         tic_toc.Stop();
+         cout << "Real time spent reordering: " <<
+              tic_toc.RealTime() << "\n";
+      }
+   }
+
+   // 6. Determine the list of true (i.e. conforming) essential boundary dofs.
+   //    In this example, the boundary conditions are defined by marking all
+   //    the boundary attributes from the mesh as essential (Dirichlet) and
+   //    converting them to a list of true dofs.
+   Array<int> ess_tdof_list;
+   if (mesh->bdr_attributes.Size())
+   {
+      Array<int> ess_bdr(mesh->bdr_attributes.Max());
+      ess_bdr = 1;
+      fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+
+   Array<int> ess_pc_tdof_list(ess_tdof_list.Size());
+
+   if (permute == 1)
+   {
+      for (int i = 0; i < ess_tdof_list.Size(); i++)
+      {
+         ess_pc_tdof_list.operator[](i) = inv_reordering->operator[](
+                                             ess_tdof_list.operator[](i));
+      }
+   }
+   else
+   {
+
+      Array<int> ess_bdr(mesh->bdr_attributes.Max());
+      ess_bdr = 1;
+      fespace->GetEssentialTrueDofs(ess_bdr, ess_pc_tdof_list);
+   }
+
+   // Array<int> ess_bdr(mesh->bdr_attributes.Max());
+   // ess_bdr = 1;
+   // fespace->GetEssentialTrueDofs(ess_bdr, ess_pc_tdof_list);
+
+   // 7. Set up the linear form b(.) which corresponds to the right-hand side
+   // of
+   //    the FEM linear system, which in this case is (1,phi_i) where phi_i are
+   //    the basis functions in the finite element fespace.
+   LinearForm *b = new LinearForm(fespace);
+
+   ConstantCoefficient one(1.0);
+   b->AddDomainIntegrator(new DomainLFIntegrator(one));
+   b->Assemble();
+
+   // 8. Define the solution vector x as a finite element grid function
+   //    corresponding to fespace. Initialize x with initial guess of zero,
+   //    which satisfies the boundary conditions.
+   GridFunction x(fespace);
+   x = 0.0;
+
+   // 9. Set up the bilinear form a(.,.) on the finite element space
+   //    corresponding to the Laplacian operator -Delta, by adding the
+   //    Diffusion domain integrator.
+   BilinearForm *a = new BilinearForm(fespace);
+   BilinearForm *a_pc = NULL;
+   if (pc) { a_pc = new BilinearForm(fespace_lor); }
+   if (pa)
+   {
+      a->SetAssemblyLevel(AssemblyLevel::PARTIAL);
+   }
+   a->AddDomainIntegrator(new DiffusionIntegrator(one));
+
+   // 10. Assemble the bilinear form and the corresponding linear system,
+   //     applying any necessary transformations such as: eliminating boundary
+   //     conditions, applying conforming constraints for non-conforming AMR,
+   //     static condensation, etc.
+   if (static_cond)
+   {
+      a->EnableStaticCondensation();
+   }
+   a->Assemble();
+
+   OperatorPtr A;
+   Vector B, X;
+
+   a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);
+
+   // 11. Solve the linear system A X = B.
+   double it_time = 0.;
+   int total_its = 0;
+
+   SparseMatrix A_pc;
+   if (pc)
+   {
+
+      tic_toc.Clear();
+      tic_toc.Start();
+
+      a_pc->AddDomainIntegrator(new DiffusionIntegrator(one));
+      a_pc->UsePrecomputedSparsity();
+      a_pc->Assemble();
+
+      tic_toc.Stop();
+      cout << "Real time assembling A_pc SparseMatrix: " <<
+           tic_toc.RealTime() << "\n";
+
+      tic_toc.Clear();
+      tic_toc.Start();
+
+      a_pc->FormSystemMatrix(ess_pc_tdof_list, A_pc);
+
+      tic_toc.Stop();
+      cout << "Real time adjusting A_pc for essential BC: " <<
+           tic_toc.RealTime() << "\n";
+
+
+      if (permute == 2)
+      {
+         tic_toc.Clear();
+         tic_toc.Start();
+
+         Array<int> perm(fespace_lor->GetTrueVSize());
+         SparseMatrix A_pc_tmp(A_pc);
+         MinimumDiscardedFillOrdering(A_pc_tmp, perm);
+         inv_reordering = new Array<int>(fespace_lor->GetTrueVSize());
+         for (int i=0; i<perm.Size(); ++i)
+         {
+            (*inv_reordering)[perm[i]] = i;
+         }
+         PermuteSparseMatrix(A_pc, *inv_reordering);
+
+         tic_toc.Stop();
+         cout << "Real time spent reordering: " <<
+              tic_toc.RealTime() << "\n";
+      }
+
+
+      if (pc_choice == GKO_BLOCK_JACOBI)
+      {
+
+         // Create Ginkgo Jacobi preconditioner
+         if (permute)
+         {
+            tic_toc.Clear();
+            tic_toc.Start();
+            GinkgoWrappers::GinkgoJacobiPreconditioner M(executor, A_pc, *inv_reordering,
+                                                         pc_storage_opt,
+                                                         pc_acc, pc_max_bs);
+            tic_toc.Stop();
+            cout << "Real time creating Ginkgo BlockJacobi preconditioner: " <<
+                 tic_toc.RealTime() << "\n";
+
+            // Use preconditioned CG
+            total_its = pcg_solve(*A, M, B, X, 0, X.Size(), 1e-12, 0.0, it_time);
+
+            cout << "Real time in PCG: " << it_time << "\n";
+         }
+         else
+         {
+            tic_toc.Clear();
+            tic_toc.Start();
+            GinkgoWrappers::GinkgoJacobiPreconditioner M(executor, A_pc, pc_storage_opt,
+                                                         pc_acc, pc_max_bs);
+            tic_toc.Stop();
+            cout << "Real time creating Ginkgo BlockJacobi preconditioner: " <<
+                 tic_toc.RealTime() << "\n";
+
+            // Use preconditioned CG
+            total_its = pcg_solve(*A, M, B, X, 0, X.Size(), 1e-12, 0.0, it_time);
+
+            cout << "Real time in PCG: " << it_time << "\n";
+         }
+
+      }
+      else if (pc_choice == GKO_ILU || pc_choice == GKO_ILU_ISAI)
+      {
+
+         // Create Ginkgo ILU preconditioner
+
+         if (permute)
+         {
+
+            tic_toc.Clear();
+            tic_toc.Start();
+
+            GinkgoWrappers::GinkgoIluPreconditioner M(executor, A_pc, *inv_reordering,
+                                                      trisolve_type, isai_sparsity_power, par_ilu_its);
+
+            tic_toc.Stop();
+            cout << "Real time creating Ginkgo Ilu preconditioner: " <<
+                 tic_toc.RealTime() << "\n";
+
+            // Use preconditioned CG
+            total_its = pcg_solve(*A, M, B, X, 0, X.Size(), 1e-12, 0.0, it_time);
+
+            cout << "Real time in PCG: " << it_time << "\n";
+
+         }
+         else
+         {
+
+            tic_toc.Clear();
+            tic_toc.Start();
+
+            GinkgoWrappers::GinkgoIluPreconditioner M(executor, A_pc, trisolve_type,
+                                                      isai_sparsity_power, par_ilu_its);
+
+            tic_toc.Stop();
+            cout << "Real time creating Ginkgo Ilu preconditioner: " <<
+                 tic_toc.RealTime() << "\n";
+
+            // Use preconditioned CG
+            total_its = pcg_solve(*A, M, B, X, 0, X.Size(), 1e-12, 0.0, it_time);
+
+            cout << "Real time in PCG: " << it_time << "\n";
+
+         }
+      }
+      else if (pc_choice == GKO_CUILU || pc_choice == GKO_CUILU_ISAI)
+      {
+
+         // Create Ginkgo CuILU preconditioner (uses cuSPARSE for factorization)
+
+         if (permute)
+         {
+
+            tic_toc.Clear();
+            tic_toc.Start();
+
+            GinkgoWrappers::GinkgoCuIluPreconditioner M(executor, A_pc, *inv_reordering,
+                                                      trisolve_type, isai_sparsity_power);
+
+            tic_toc.Stop();
+            cout << "Real time creating Ginkgo CuIlu preconditioner: " <<
+                 tic_toc.RealTime() << "\n";
+
+            // Use preconditioned CG
+            total_its = pcg_solve(*A, M, B, X, 0, X.Size(), 1e-12, 0.0, it_time);
+
+            cout << "Real time in PCG: " << it_time << "\n";
+
+         }
+         else
+         {
+
+            tic_toc.Clear();
+            tic_toc.Start();
+
+            GinkgoWrappers::GinkgoCuIluPreconditioner M(executor, A_pc, trisolve_type,
+                                                      isai_sparsity_power);
+
+            tic_toc.Stop();
+            cout << "Real time creating Ginkgo CuIlu preconditioner: " <<
+                 tic_toc.RealTime() << "\n";
+
+            // Use preconditioned CG
+            total_its = pcg_solve(*A, M, B, X, 0, X.Size(), 1e-12, 0.0, it_time);
+
+            cout << "Real time in PCG: " << it_time << "\n";
+
+         }
+      }
+      else if (pc_choice == MFEM_GS)
+      {
+
+         // Create MFEM preconditioner
+         tic_toc.Clear();
+         tic_toc.Start();
+
+         GSSmoother M(A_pc);
+
+         tic_toc.Stop();
+         cout << "Real time creating MFEM GS preconditioner: " <<
+              tic_toc.RealTime() << "\n";
+
+         // Use preconditioned CG
+         total_its = pcg_solve(*A, M, B, X, 0, X.Size(), 1e-12, 0.0, it_time);
+
+         cout << "Real time in PCG: " << it_time << "\n";
+
+      }
+      else if (pc_choice == MFEM_UMFPACK)
+      {
+
+#ifdef MFEM_USE_SUITESPARSE
+         // Create MFEM preconditioner
+         tic_toc.Clear();
+         tic_toc.Start();
+
+         UMFPackSolver M;
+         M.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
+         M.SetOperator(A_pc);
+
+         tic_toc.Stop();
+         cout << "Real time creating MFEM UMFPACK preconditioner: " <<
+              tic_toc.RealTime() << "\n";
+
+         // Use preconditioned CG
+         total_its = pcg_solve(*A, M, B, X, 0, X.Size(), 1e-12, 0.0, it_time);
+
+         cout << "Real time in PCG: " << it_time << "\n";
+#endif
+      }
+   }
+   else
+   {
+
+      total_its = cg_solve(*A, B, X, 0, X.Size(), 1e-12, 0.0, it_time);
+
+      cout << "Real time in CG: " << it_time << "\n";
+   }
+
+   cout << "Total iterations: " << total_its << "\n";
+   cout << "Avg time per iteration: " << it_time/double(total_its) << "\n";
+
+   // 12. Recover the solution as a finite element grid function.
+   a->RecoverFEMSolution(X, *b, x);
+
+   // 13. Save the refined mesh and the solution. This output can be viewed
+   // later
+   //     using GLVis: "glvis -m refined.mesh -g sol.gf".
+
+   if (output_sol)
+   {
+      ofstream mesh_ofs("refined.mesh");
+      mesh_ofs.precision(8);
+      mesh->Print(mesh_ofs);
+      ofstream sol_ofs("sol.gf");
+      sol_ofs.precision(8);
+      x.Save(sol_ofs);
+   }
+
+   if (pc && output_pc)
+   {
+      ofstream mesh_lor_ofs("lor-refined.mesh");
+      mesh_lor_ofs.precision(8);
+      mesh_lor->Print(mesh_lor_ofs);
+
+      ofstream apc_lor_ofs("lor-mat.dat");
+      mesh_lor_ofs.precision(8);
+      A_pc.PrintCSR(apc_lor_ofs);
+   }
+
+   // 14. Send the solution by socket to a GLVis server.
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int visport = 19916;
+      socketstream sol_sock(vishost, visport);
+      sol_sock.precision(8);
+      sol_sock << "solution\n" << *mesh << x << flush;
+   }
+
+   // 15. Free the used memory.
+   delete a;
+   delete b;
+   delete fespace;
+   delete fespace_lor;
+   delete fec_lor;
+   delete mesh_lor;
+   if (order > 0)
+   {
+      delete fec;
+   }
+   delete mesh;
+}
+

examples/ginkgo/multigridpc.hpp

@@ -0,0 +1,217 @@
+#ifndef __MULTIGRIDPC_HPP__
+#define __MULTIGRIDPC_HPP__
+
+#include "mfem.hpp"
+#include <memory>
+
+namespace mfem
+{
+
+struct SolverConfig
+{
+   enum SolverType
+   {
+     JACOBI = 0,
+     AMGX = 1,
+     CHEBYSHEV = 2,
+     GINKGO_CUIC = 3,
+     GINKGO_CUIC_ISAI = 4 
+   };
+   SolverType type;
+   SolverType smoother_type;
+   AssemblyLevel upper_level_asm;
+   const char *amgx_config;
+   std::shared_ptr<gko::Executor> gko_exec;
+#ifdef MFEM_SIMPLEX_LOR
+   bool simplex_lor = true;
+#endif
+   SolverConfig(SolverType type_, SolverType sm_type_, AssemblyLevel upper_asm_, 
+                const char *amgx_config_, std::shared_ptr<gko::Executor> gko_exec_) : type(type_),
+     smoother_type(sm_type_),
+     upper_level_asm(upper_asm_) 
+     {
+       amgx_config = amgx_config_;
+       gko_exec = gko_exec_; 
+     } 
+};
+
+struct MGRefinement
+{  
+   enum Type { P_MG, H_MG };
+   Type type;
+   int order;
+   MGRefinement(Type type_, int order_) : type(type_), order(order_) { }
+   static MGRefinement p(int order_) { return MGRefinement(P_MG, order_); }
+   static MGRefinement h() { return MGRefinement(H_MG, 0); }
+};
+
+bool NeedsLOR(SolverConfig config)
+{
+   switch (config.type)
+   {
+      case SolverConfig::GINKGO_CUIC:
+      case SolverConfig::GINKGO_CUIC_ISAI:
+         return true;
+      default:
+         return false;
+   }
+}
+
+struct DiffusionMultigrid : GeometricMultigrid
+{
+   Coefficient &coeff;
+   OperatorPtr A_coarse;
+
+   DiffusionMultigrid(
+      FiniteElementSpaceHierarchy& hierarchy,
+      Coefficient &coeff_,
+      Array<int>& ess_bdr,
+      SolverConfig solver_config)
+       : GeometricMultigrid(hierarchy), coeff(coeff_)
+    {
+       ConstructCoarseOperatorAndSolver(
+          solver_config, hierarchy.GetFESpaceAtLevel(0), ess_bdr);
+       int nlevels = hierarchy.GetNumLevels();
+       for (int i=1; i<nlevels; ++i)
+       {
+          ConstructOperatorAndSmoother(solver_config, hierarchy.GetFESpaceAtLevel(i), ess_bdr);
+       }
+    }
+
+   void ConstructBilinearForm(
+      FiniteElementSpace &fespace, Array<int> &ess_bdr, AssemblyLevel asm_lvl)
+    {
+       BilinearForm* form = new BilinearForm(&fespace);
+       form->SetAssemblyLevel(asm_lvl);
+       form->SetDiagonalPolicy(DIAG_ONE);
+       form->AddDomainIntegrator(new DiffusionIntegrator(coeff));
+       form->Assemble();
+       bfs.Append(form);
+    
+       essentialTrueDofs.Append(new Array<int>());
+       fespace.GetEssentialTrueDofs(ess_bdr, *essentialTrueDofs.Last());
+    }
+
+   void ConstructOperatorAndSmoother(SolverConfig solver_config,
+      FiniteElementSpace& fespace, Array<int>& ess_bdr)
+    {
+       ConstructBilinearForm(fespace, ess_bdr, solver_config.upper_level_asm);
+    
+       OperatorPtr opr;
+       bfs.Last()->FormSystemMatrix(*essentialTrueDofs.Last(), opr);
+       opr.SetOperatorOwner(false);
+   
+       switch (solver_config.smoother_type)
+       {
+         case SolverConfig::CHEBYSHEV: 
+         {
+           Vector diag(fespace.GetTrueVSize());
+           bfs.Last()->AssembleDiagonal(diag);
+    
+           Solver* smoother = new OperatorChebyshevSmoother(
+              opr.Ptr(), diag, *essentialTrueDofs.Last(), 2);
+    
+           if (solver_config.upper_level_asm == AssemblyLevel::PARTIAL) 
+           {
+             AddLevel(opr.Ptr(), smoother, true, true);
+           }
+           else
+           {
+             AddLevel(opr.Ptr(), smoother, false, true);
+           }
+           break;
+         }
+         case SolverConfig::GINKGO_CUIC:
+         {
+            SparseMatrix *A_lvl = dynamic_cast<SparseMatrix*>(opr.Ptr());
+            Solver *smoother = new GinkgoWrappers::GinkgoCuIcPreconditioner(
+                                              solver_config.gko_exec, *A_lvl,
+                                                        "exact", 1);
+           if (solver_config.upper_level_asm == AssemblyLevel::PARTIAL) 
+           {
+             AddLevel(opr.Ptr(), smoother, true, true);
+           }
+           else
+           {
+             AddLevel(opr.Ptr(), smoother, false, true);
+           }
+           break;
+         }
+         case SolverConfig::GINKGO_CUIC_ISAI:
+         {
+            SparseMatrix *A_lvl = dynamic_cast<SparseMatrix*>(opr.Ptr());
+            Solver *smoother = new GinkgoWrappers::GinkgoCuIcPreconditioner(
+                                              solver_config.gko_exec, *A_lvl,
+                                                        "isai", 1);
+           if (solver_config.upper_level_asm == AssemblyLevel::PARTIAL) 
+           {
+             AddLevel(opr.Ptr(), smoother, true, true);
+           }
+           else
+           {
+             AddLevel(opr.Ptr(), smoother, false, true);
+           }
+           break;
+         }
+       }
+      
+    }
+
+   void ConstructCoarseOperatorAndSolver(
+      SolverConfig config, FiniteElementSpace& fespace, Array<int>& ess_bdr)
+    {
+       ConstructBilinearForm(fespace, ess_bdr, AssemblyLevel::LEGACYFULL);
+       BilinearForm &a = *bfs.Last();
+       Array<int> &ess_dofs = *essentialTrueDofs.Last();
+    
+       bfs.Last()->FormSystemMatrix(*essentialTrueDofs.Last(), A_coarse);
+    
+       OperatorPtr A_prec;
+       {
+          A_prec = A_coarse;
+       }
+    
+       Solver *coarse_solver;
+       switch (config.type)
+       {
+          case SolverConfig::JACOBI:
+             coarse_solver = new OperatorJacobiSmoother(a, ess_dofs);
+             break;
+    #ifdef MFEM_USE_AMGX
+          case SolverConfig::AMGX:
+          {
+             AmgXSolver *amg = new AmgXSolver;
+             amg->ReadParameters(config.amgx_config, AmgXSolver::EXTERNAL);
+             amg->InitSerial();
+             amg->SetOperator(*A_prec.As<SparseMatrix>());
+             coarse_solver = amg;
+             break;
+          }
+    #endif
+          case SolverConfig::GINKGO_CUIC:
+          {
+            SparseMatrix *A_lvl = dynamic_cast<SparseMatrix*>(A_prec.Ptr());
+            Solver *gko_solver = new GinkgoWrappers::GinkgoCuIcPreconditioner(
+                                              config.gko_exec, *A_lvl,
+                                                        "exact", 1);
+            coarse_solver = gko_solver;
+            break;
+          }
+          case SolverConfig::GINKGO_CUIC_ISAI:
+          {
+            SparseMatrix *A_lvl = dynamic_cast<SparseMatrix*>(A_prec.Ptr());
+            Solver *gko_solver = new GinkgoWrappers::GinkgoCuIcPreconditioner(
+                                              config.gko_exec, *A_lvl,
+                                                        "isai", 1);
+            coarse_solver = gko_solver;
+            break;
+          }
+          default:
+             MFEM_ABORT("Not available.")
+       }
+    
+       AddLevel(A_coarse.Ptr(), coarse_solver, false, true);
+    }
+};
+}
+#endif                                                             

fem/bilinearform.hpp

@@ -299,8 +299,17 @@ public:
    /// Returns a reference to the sparse matrix:  \f$ M \f$
    SparseMatrix &SpMat()
    {
-      MFEM_VERIFY(mat, "mat is NULL and can't be dereferenced");
-      return *mat;
+      FABilinearFormExtension *fa_ext = dynamic_cast<FABilinearFormExtension*>
+                                        (this->ext);
+      if (fa_ext)
+      {
+         return fa_ext->SpMat();
+      }
+      else
+      {
+         MFEM_VERIFY(mat, "mat is NULL and can't be dereferenced");
+         return *mat;
+      }
    }
 
    /**  @brief Nullifies the internal matrix \f$ M \f$ and returns a pointer

fem/bilinearform_ext.hpp

@@ -129,6 +129,10 @@ public:
    void Assemble();
    void Mult(const Vector &x, Vector &y) const;
    void MultTranspose(const Vector &x, Vector &y) const;
+   SparseMatrix &SpMat()
+   {
+      return mat;
+   }
 };
 
 /// Data and methods for matrix-free bilinear forms

fem/fespace.cpp

@@ -2806,9 +2806,9 @@ L2ProjectionGridTransfer::L2Projection::L2Projection(
    Vector shape_ho(ndof_ho);
    Vector shape_lor(ndof_lor);
 
-   const Geometry::Type geom = fe_ho->GetGeomType();
-   const DenseTensor &pmats = cf_tr.point_matrices[geom];
-   emb_tr.SetIdentityTransformation(geom);
+   const Geometry::Type lor_geom = fe_lor->GetGeomType();
+   const DenseTensor &pmats = cf_tr.point_matrices[lor_geom];
+   emb_tr.SetIdentityTransformation(lor_geom);
 
    for (int iho=0; iho<nel_ho; ++iho)
    {
@@ -2833,7 +2833,7 @@ L2ProjectionGridTransfer::L2Projection::L2Projection(
          emb_tr.SetPointMat(pmats(cf_tr.embeddings[ilor].matrix));
 
          int order = fe_lor->GetOrder() + fe_ho->GetOrder() + el_tr->OrderW();
-         const IntegrationRule *ir = &IntRules.Get(geom, order);
+         const IntegrationRule *ir = &IntRules.Get(lor_geom, order);
          M_mixed_el = 0.0;
          for (int i = 0; i < ir->GetNPoints(); i++)
          {

linalg/ginkgo.cpp

@@ -276,8 +276,8 @@ CGSolver::CGSolver(
                                ATOLERANCE)
 {
    using cg = gko::solver::Cg<double>;
-   this->solver_gen =
-      cg::build().with_criteria(this->combined_factory).on(this->executor);
+   //   this->solver_gen =
+   //      cg::build().with_criteria(this->combined_factory).on(this->executor);
 }
 
 CGSolver::CGSolver(
@@ -292,10 +292,10 @@ CGSolver::CGSolver(
                                ATOLERANCE)
 {
    using cg         = gko::solver::Cg<double>;
-   this->solver_gen = cg::build()
-                      .with_criteria(this->combined_factory)
-                      .with_preconditioner(preconditioner)
-                      .on(this->executor);
+   //   this->solver_gen = cg::build()
+   //                      .with_criteria(this->combined_factory)
+   //                      .with_preconditioner(preconditioner)
+   //                      .on(this->executor);
 }
 
 
@@ -311,9 +311,9 @@ BICGSTABSolver::BICGSTABSolver(
                                ATOLERANCE)
 {
    using bicgstab   = gko::solver::Bicgstab<double>;
-   this->solver_gen = bicgstab::build()
-                      .with_criteria(this->combined_factory)
-                      .on(this->executor);
+   //   this->solver_gen = bicgstab::build()
+   //                      .with_criteria(this->combined_factory)
+   //                      .on(this->executor);
 }
 
 BICGSTABSolver::BICGSTABSolver(
@@ -328,10 +328,10 @@ BICGSTABSolver::BICGSTABSolver(
                                ATOLERANCE)
 {
    using bicgstab   = gko::solver::Bicgstab<double>;
-   this->solver_gen = bicgstab::build()
-                      .with_criteria(this->combined_factory)
-                      .with_preconditioner(preconditioner)
-                      .on(this->executor);
+   //   this->solver_gen = bicgstab::build()
+   //                      .with_criteria(this->combined_factory)
+   //                      .with_preconditioner(preconditioner)
+   //                      .on(this->executor);
 }
 
 
@@ -347,8 +347,8 @@ CGSSolver::CGSSolver(
                                ATOLERANCE)
 {
    using cgs = gko::solver::Cgs<double>;
-   this->solver_gen =
-      cgs::build().with_criteria(this->combined_factory).on(this->executor);
+   //   this->solver_gen =
+   //      cgs::build().with_criteria(this->combined_factory).on(this->executor);
 }
 
 CGSSolver::CGSSolver(
@@ -363,10 +363,10 @@ CGSSolver::CGSSolver(
                                ATOLERANCE)
 {
    using cgs        = gko::solver::Cgs<double>;
-   this->solver_gen = cgs::build()
-                      .with_criteria(this->combined_factory)
-                      .with_preconditioner(preconditioner)
-                      .on(this->executor);
+   //   this->solver_gen = cgs::build()
+   //                      .with_criteria(this->combined_factory)
+   //                      .with_preconditioner(preconditioner)
+   //                      .on(this->executor);
 }
 
 
@@ -382,8 +382,8 @@ FCGSolver::FCGSolver(
                                ATOLERANCE)
 {
    using fcg = gko::solver::Fcg<double>;
-   this->solver_gen =
-      fcg::build().with_criteria(this->combined_factory).on(this->executor);
+   //   this->solver_gen =
+   //      fcg::build().with_criteria(this->combined_factory).on(this->executor);
 }
 
 FCGSolver::FCGSolver(
@@ -398,10 +398,10 @@ FCGSolver::FCGSolver(
                                ATOLERANCE)
 {
    using fcg        = gko::solver::Fcg<double>;
-   this->solver_gen = fcg::build()
-                      .with_criteria(this->combined_factory)
-                      .with_preconditioner(preconditioner)
-                      .on(this->executor);
+   //   this->solver_gen = fcg::build()
+   //                      .with_criteria(this->combined_factory)
+   //                      .with_preconditioner(preconditioner)
+   //                      .on(this->executor);
 }
 
 
@@ -417,10 +417,10 @@ GMRESSolver::GMRESSolver(
                                ATOLERANCE)
 {
    using gmres      = gko::solver::Gmres<double>;
-   this->solver_gen = gmres::build()
-                      .with_krylov_dim(m)
-                      .with_criteria(this->combined_factory)
-                      .on(this->executor);
+   //   this->solver_gen = gmres::build()
+   //                      .with_krylov_dim(m)
+   //                      .with_criteria(this->combined_factory)
+   //                      .on(this->executor);
 }
 
 GMRESSolver::GMRESSolver(
@@ -435,11 +435,11 @@ GMRESSolver::GMRESSolver(
                                ATOLERANCE)
 {
    using gmres      = gko::solver::Gmres<double>;
-   this->solver_gen = gmres::build()
-                      .with_krylov_dim(m)
-                      .with_criteria(this->combined_factory)
-                      .with_preconditioner(preconditioner)
-                      .on(this->executor);
+   //   this->solver_gen = gmres::build()
+   //                      .with_krylov_dim(m)
+   //                      .with_criteria(this->combined_factory)
+   //                      .with_preconditioner(preconditioner)
+   //                      .on(this->executor);
 }
 
 
@@ -455,8 +455,8 @@ IRSolver::IRSolver(
                                ATOLERANCE)
 {
    using ir = gko::solver::Ir<double>;
-   this->solver_gen =
-      ir::build().with_criteria(this->combined_factory).on(this->executor);
+   //   this->solver_gen =
+   //      ir::build().with_criteria(this->combined_factory).on(this->executor);
 }
 
 IRSolver::IRSolver(
@@ -471,10 +471,105 @@ IRSolver::IRSolver(
                                ATOLERANCE)
 {
    using ir         = gko::solver::Ir<double>;
-   this->solver_gen = ir::build()
-                      .with_criteria(this->combined_factory)
-                      .with_solver(inner_solver)
-                      .on(this->executor);
+   //   this->solver_gen = ir::build()
+   //                      .with_criteria(this->combined_factory)
+   //                      .with_solver(inner_solver)
+   //                      .on(this->executor);
+}
+
+
+/* ---------------------- GinkgoPreconditioner ------------------------ */
+void GinkgoPreconditionerBase::Mult(const Vector &x, Vector &y) const
+{
+
+   if (!iterative_mode)
+   {
+      y = 0.0;
+   }
+
+   //Create Ginkgo wrapped-vectors
+   bool on_device = false;
+   if (exec_ != exec_->get_master())
+   {
+      on_device = true;
+   }
+   using vec = gko::matrix::Dense<double>;
+   auto gko_x = vec::create(exec_, gko::dim<2> {x.Size(), 1},
+                            gko::Array<double>::view(exec_,
+                                                     x.Size(), const_cast<double *>(
+                                                        x.Read(on_device))), 1);
+   auto gko_y = vec::create(exec_, gko::dim<2> {y.Size(), 1},
+                            gko::Array<double>::view(exec_,
+                                                     y.Size(), y.ReadWrite(on_device)), 1);
+
+   if (permute_)
+   {
+
+      auto gko_x_perm = vec::create(exec_, gko::dim<2> {x.Size(), 1});
+      auto gko_y_perm = vec::create(exec_, gko::dim<2> {y.Size(), 1});
+
+      vec_permute_->apply(gko::lend(gko_x), gko::lend(gko_x_perm));
+      vec_permute_->apply(gko::lend(gko_y), gko::lend(gko_y_perm));
+
+      gko_precond_.get()->apply(gko::lend(gko_x_perm), gko::lend(gko_y_perm));
+
+      vec_inv_permute_->apply(gko::lend(gko_y_perm), gko::lend(gko_y));
+
+      // Reset MFEM Vector to use Ginkgo output ?
+      //   y = gko_y->get_values();
+
+   }
+   else
+   {
+
+      gko_precond_.get()->apply(gko::lend(gko_x), gko::lend(gko_y));
+
+   }
+}
+
+// TEMP!  This only works for symmetric operators
+void GinkgoPreconditionerBase::MultTranspose(const Vector &x, Vector &y) const
+{
+  this->Mult(x, y);
+}
+
+void GinkgoPreconditionerBase::SetOperator(const Operator &op)
+{
+
+   // Only accept SparseMatrix for this type (see SparseSmoother::SetOperator)
+   SparseMatrix *op_mat = const_cast<SparseMatrix*>(
+                             dynamic_cast<const SparseMatrix*>(&op));
+   if (op_mat == NULL)
+   {
+      mfem_error("GinkgoPreconditioner::SetOperator : not a SparseMatrix!");
+   }
+   height = op_mat->Height();
+   width = op_mat->Width();
+
+   // Release current preconditioner
+   gko_precond_.release();
+
+   bool on_device = false;
+   if (exec_ != exec_->get_master())
+   {
+      on_device = true;
+   }
+
+   using mtx = gko::matrix::Csr<double, int>;
+   auto gko_sparse = mtx::create(
+                        exec_, gko::dim<2>(op_mat->Height(), op_mat->Width()),
+                        gko::Array<double>::view(exec_,
+                                                 op_mat->NumNonZeroElems(),
+                                                 op_mat->ReadWriteData(on_device)),
+                        gko::Array<int>::view(exec_,
+                                              op_mat->NumNonZeroElems(),
+                                              op_mat->ReadWriteJ(on_device)),
+                        gko::Array<int>::view(exec_,
+                                              op_mat->Height() + 1,
+                                              op_mat->ReadWriteI(on_device)));
+   gko_precond_ = gko_precond_factory_.get()->generate(
+                     gko::give(gko_sparse));
+
 }
 
 

linalg/ginkgo.hpp

@@ -640,6 +640,725 @@ public:
 };
 
 
+class GinkgoPreconditionerBase : public Solver
+{
+protected:
+   GinkgoPreconditionerBase(std::shared_ptr<const gko::Executor> exec,
+                            SparseMatrix &a,
+                            bool iter_mode=false)
+      : Solver(a.Height(), a.Width(), iter_mode)
+   {
+      exec_ = std::move(exec);
+      permute_ = false;
+   }
+
+   GinkgoPreconditionerBase(std::shared_ptr<const gko::Executor> exec,
+                            SparseMatrix &a, Array<int> &inv_permutation_indices,
+                            bool iter_mode=false)
+      : Solver(a.Height(), a.Width(), iter_mode)
+   {
+      exec_ = std::move(exec);
+
+      bool on_device = false;
+      if (exec_->get_master() != exec_)
+      {
+         on_device = true;
+      }
+      permute_ = true;
+      auto gko_inv_perm_ind = gko::Array<int>::view(
+                                 exec_,
+                                 inv_permutation_indices.Size(),
+                                 inv_permutation_indices.ReadWrite(
+                                    on_device));
+      // Note the "forward" permutation uses the inverse flag because
+      //  the indices are for the inverse permutation as defined by Ginkgo
+      vec_permute_ = gko::matrix::Permutation<int>::create(
+                        exec_, gko::dim<2> {inv_permutation_indices.Size()},
+                        gko_inv_perm_ind,
+                        gko::matrix::row_permute | gko::matrix::inverse_permute);
+      vec_inv_permute_ = gko::matrix::Permutation<int>::create(
+                            exec_, gko::dim<2> {inv_permutation_indices.Size()},
+                            gko_inv_perm_ind,
+                            gko::matrix::row_permute);
+   }
+
+public:
+   std::shared_ptr<const gko::Executor> get_exec() {return this->exec_; }
+   const gko::LinOpFactory* get_gko_precond_factory()
+   {
+      return this->gko_precond_factory_.get();
+   }
+   gko::LinOp* get_gko_precond() {return this->gko_precond_.get(); }
+
+   virtual void Mult(const Vector &x, Vector &y) const;
+   virtual void MultTranspose(const Vector &x, Vector &y) const;
+   virtual void SetOperator(const Operator &op);
+
+protected:
+   std::shared_ptr<const gko::Executor> exec_;
+   std::unique_ptr<const gko::LinOpFactory> gko_precond_factory_;
+   std::unique_ptr<gko::LinOp> gko_precond_;
+   bool permute_;
+   std::unique_ptr<gko::matrix::Permutation<int>> vec_permute_;
+   std::unique_ptr<gko::matrix::Permutation<int>> vec_inv_permute_;
+
+};
+
+
+class GinkgoJacobiPreconditioner : public GinkgoPreconditionerBase
+{
+public:
+   GinkgoJacobiPreconditioner(std::shared_ptr<const gko::Executor> exec,
+                              SparseMatrix &a,
+                              const char *storage_opt="none",
+                              const double accuracy=1.e-1,
+                              const int max_block_size=32,
+                              bool iter_mode=false)
+      : GinkgoPreconditionerBase(exec, a, iter_mode)
+   {
+
+      bool on_device = false;
+      if (exec->get_master() != exec)
+      {
+         on_device = true;
+      }
+
+      using mtx = gko::matrix::Csr<double, int>;
+      const int nnz =  a.GetMemoryData().Capacity();
+      auto gko_sparse = mtx::create(
+                           exec, gko::dim<2>(a.Height(), a.Width()),
+                           gko::Array<double>::view(exec,
+                                                    nnz,
+                                                    a.ReadWriteData(on_device)),
+                           gko::Array<int>::view(exec,
+                                                 nnz,
+                                                 a.ReadWriteJ(on_device)),
+                           gko::Array<int>::view(exec, a.Height() + 1,
+                                                 a.ReadWriteI(on_device)));
+
+      gko_sparse->sort_by_column_index();
+
+      if (storage_opt == "auto")
+      {
+         gko_precond_factory_ = gko::preconditioner::Jacobi<double, int>::build()
+                                .with_storage_optimization(
+                                   gko::precision_reduction::autodetect())
+                                .with_accuracy(accuracy)
+                                .with_max_block_size(max_block_size)
+                                .on(exec);
+      }
+      else
+      {
+         gko_precond_factory_ = gko::preconditioner::Jacobi<double, int>::build()
+                                .with_storage_optimization(
+                                   gko::precision_reduction(0, 0))
+                                .with_accuracy(accuracy)
+                                .with_max_block_size(max_block_size)
+                                .on(exec);
+      }
+
+      gko_precond_ = gko_precond_factory_.get()->generate(
+                        gko::give(gko_sparse));
+   }
+
+   GinkgoJacobiPreconditioner(std::shared_ptr<const gko::Executor> exec,
+                              SparseMatrix &a,
+                              Array<int> &inv_permutation_indices,
+                              const char *storage_opt="none",
+                              const double accuracy=1.e-1,
+                              const int max_block_size=32,
+                              bool iter_mode=false)
+      : GinkgoPreconditionerBase(exec, a, inv_permutation_indices, iter_mode)
+   {
+
+      bool on_device = false;
+      if (exec->get_master() != exec)
+      {
+         on_device = true;
+      }
+
+      using mtx = gko::matrix::Csr<double, int>;
+      const int nnz =  a.GetMemoryData().Capacity();
+      auto gko_sparse = mtx::create(
+                           exec, gko::dim<2>(a.Height(), a.Width()),
+                           gko::Array<double>::view(exec,
+                                                    nnz,
+                                                    a.ReadWriteData(on_device)),
+                           gko::Array<int>::view(exec,
+                                                 nnz,
+                                                 a.ReadWriteJ(on_device)),
+                           gko::Array<int>::view(exec, a.Height() + 1,
+                                                 a.ReadWriteI(on_device)));
+
+
+      gko_sparse->sort_by_column_index();
+
+      if (storage_opt == "auto")
+      {
+         gko_precond_factory_ = gko::preconditioner::Jacobi<double, int>::build()
+                                .with_storage_optimization(
+                                   gko::precision_reduction::autodetect())
+                                .with_accuracy(accuracy)
+                                .with_max_block_size(max_block_size)
+                                .on(exec);
+      }
+      else
+      {
+         gko_precond_factory_ = gko::preconditioner::Jacobi<double, int>::build()
+                                .with_storage_optimization(
+                                   gko::precision_reduction(0, 0))
+                                .with_accuracy(accuracy)
+                                .with_max_block_size(max_block_size)
+                                .on(exec);
+      }
+
+      gko_precond_ = gko_precond_factory_.get()->generate(
+                        gko::give(gko_sparse));
+   }
+};
+
+class GinkgoIcPreconditioner : public GinkgoPreconditionerBase
+{
+public:
+   GinkgoIcPreconditioner(std::shared_ptr<const gko::Executor> exec,
+                           SparseMatrix &a, const char *trisolve_type = "exact",
+                           int sparsity_power=1,
+                           int par_ic_its=0,
+                           bool iter_mode=false)
+      : GinkgoPreconditionerBase(exec, a, iter_mode)
+   {
+
+      bool on_device = false;
+      if (exec->get_master() != exec)
+      {
+         on_device = true;
+      }
+
+      using mtx = gko::matrix::Csr<double, int>;
+      const int nnz =  a.GetMemoryData().Capacity();
+      auto gko_sparse = mtx::create(
+                           exec, gko::dim<2>(a.Height(), a.Width()),
+                           gko::Array<double>::view(exec,
+                                                    nnz,
+                                                    a.ReadWriteData(on_device)),
+                           gko::Array<int>::view(exec,
+                                                 nnz,
+                                                 a.ReadWriteJ(on_device)),
+                           gko::Array<int>::view(exec, a.Height() + 1,
+                                                 a.ReadWriteI(on_device)));
+
+      using ic_fact_type = gko::factorization::ParIc<double, int>;
+      std::shared_ptr<ic_fact_type::Factory> fact_factory = std::move(
+                                                                ic_fact_type::build()
+                                                                .with_iterations(par_ic_its)
+                                                                .with_both_factors(false)
+                                                                .on(exec));
+
+
+      if (trisolve_type == "isai")
+      {
+
+         using l_solver_type = gko::preconditioner::LowerIsai<>;
+
+         std::shared_ptr<l_solver_type::Factory> l_solver_factory = std::move(
+                                                                       l_solver_type::build()
+                                                                       .with_sparsity_power(sparsity_power)
+                                                                       .on(exec));
+
+         gko_precond_factory_ = gko::preconditioner::Ic<l_solver_type>::build()
+         .with_factorization_factory(fact_factory)
+         .with_l_solver_factory(l_solver_factory)
+         .on(exec);
+
+      }
+      else
+      {
+
+         gko_precond_factory_ = gko::preconditioner::Ic<>::build()
+                                .with_factorization_factory(fact_factory)
+                                .on(exec);
+      }
+
+      gko_precond_ = gko_precond_factory_.get()->generate(
+                        gko::give(gko_sparse));
+   }
+
+   GinkgoIcPreconditioner(std::shared_ptr<const gko::Executor> exec,
+                           SparseMatrix &a, Array<int> &inv_permutation_indices,
+                           const char *trisolve_type = "exact",
+                           int sparsity_power=1,
+                           int par_ic_its=0,
+                           bool iter_mode=false)
+      : GinkgoPreconditionerBase(exec, a, inv_permutation_indices, iter_mode)
+   {
+
+      bool on_device = false;
+      if (exec->get_master() != exec)
+      {
+         on_device = true;
+      }
+
+      using mtx = gko::matrix::Csr<double, int>;
+      const int nnz =  a.GetMemoryData().Capacity();
+      auto gko_sparse = mtx::create(
+                           exec, gko::dim<2>(a.Height(), a.Width()),
+                           gko::Array<double>::view(exec,
+                                                    nnz,
+                                                    a.ReadWriteData(on_device)),
+                           gko::Array<int>::view(exec,
+                                                 nnz,
+                                                 a.ReadWriteJ(on_device)),
+                           gko::Array<int>::view(exec, a.Height() + 1,
+                                                 a.ReadWriteI(on_device)));
+
+      using ic_fact_type = gko::factorization::ParIc<double, int>;
+      std::shared_ptr<ic_fact_type::Factory> fact_factory = std::move(
+                                                                ic_fact_type::build()
+                                                                .with_iterations(par_ic_its)
+                                                                .with_both_factors(false)
+                                                                .on(exec));
+
+      if (trisolve_type == "isai")
+      {
+
+         using l_solver_type = gko::preconditioner::LowerIsai<>;
+
+         std::shared_ptr<l_solver_type::Factory> l_solver_factory = std::move(
+                                                                       l_solver_type::build()
+                                                                       .with_sparsity_power(sparsity_power)
+                                                                       .on(exec));
+
+         gko_precond_factory_ = gko::preconditioner::Ic<l_solver_type>::build()
+         .with_factorization_factory(fact_factory)
+         .with_l_solver_factory(l_solver_factory)
+         .on(exec);
+      }
+      else
+      {
+
+         gko_precond_factory_ = gko::preconditioner::Ic<>::build()
+                                .with_factorization_factory(fact_factory)
+                                .on(exec);
+      }
+
+      gko_precond_ = gko_precond_factory_.get()->generate(
+                        gko::give(gko_sparse));
+   }
+};
+
+class GinkgoIluPreconditioner : public GinkgoPreconditionerBase
+{
+public:
+   GinkgoIluPreconditioner(std::shared_ptr<const gko::Executor> exec,
+                           SparseMatrix &a, const char *trisolve_type = "exact",
+                           int sparsity_power=1,
+                           int par_ilu_its=0,
+                           bool iter_mode=false)
+      : GinkgoPreconditionerBase(exec, a, iter_mode)
+   {
+
+      bool on_device = false;
+      if (exec->get_master() != exec)
+      {
+         on_device = true;
+      }
+
+      using mtx = gko::matrix::Csr<double, int>;
+      const int nnz =  a.GetMemoryData().Capacity();
+      auto gko_sparse = mtx::create(
+                           exec, gko::dim<2>(a.Height(), a.Width()),
+                           gko::Array<double>::view(exec,
+                                                    nnz,
+                                                    a.ReadWriteData(on_device)),
+                           gko::Array<int>::view(exec,
+                                                 nnz,
+                                                 a.ReadWriteJ(on_device)),
+                           gko::Array<int>::view(exec, a.Height() + 1,
+                                                 a.ReadWriteI(on_device)));
+
+      using ilu_fact_type = gko::factorization::ParIlu<double, int>;
+      std::shared_ptr<ilu_fact_type::Factory> fact_factory = std::move(
+                                                                ilu_fact_type::build()
+                                                                .with_iterations(par_ilu_its)
+                                                                .on(exec));
+
+
+      if (trisolve_type == "isai")
+      {
+
+         using l_solver_type = gko::preconditioner::LowerIsai<>;
+         using u_solver_type = gko::preconditioner::UpperIsai<>;
+
+         std::shared_ptr<l_solver_type::Factory> l_solver_factory = std::move(
+                                                                       l_solver_type::build()
+                                                                       .with_sparsity_power(sparsity_power)
+                                                                       .on(exec));
+         std::shared_ptr<u_solver_type::Factory> u_solver_factory = std::move(
+                                                                       u_solver_type::build()
+                                                                       .with_sparsity_power(sparsity_power)
+                                                                       .on(exec));
+
+         gko_precond_factory_ = gko::preconditioner::Ilu<l_solver_type,
+         u_solver_type>::build()
+         .with_factorization_factory(fact_factory)
+         .with_l_solver_factory(l_solver_factory)
+         .with_u_solver_factory(u_solver_factory)
+         .on(exec);
+
+      }
+      else
+      {
+
+         gko_precond_factory_ = gko::preconditioner::Ilu<>::build()
+                                .with_factorization_factory(fact_factory)
+                                .on(exec);
+      }
+
+      gko_precond_ = gko_precond_factory_.get()->generate(
+                        gko::give(gko_sparse));
+   }
+
+   GinkgoIluPreconditioner(std::shared_ptr<const gko::Executor> exec,
+                           SparseMatrix &a, Array<int> &inv_permutation_indices,
+                           const char *trisolve_type = "exact",
+                           int sparsity_power=1,
+                           int par_ilu_its=0,
+                           bool iter_mode=false)
+      : GinkgoPreconditionerBase(exec, a, inv_permutation_indices, iter_mode)
+   {
+
+      bool on_device = false;
+      if (exec->get_master() != exec)
+      {
+         on_device = true;
+      }
+
+      using mtx = gko::matrix::Csr<double, int>;
+      const int nnz =  a.GetMemoryData().Capacity();
+      auto gko_sparse = mtx::create(
+                           exec, gko::dim<2>(a.Height(), a.Width()),
+                           gko::Array<double>::view(exec,
+                                                    nnz,
+                                                    a.ReadWriteData(on_device)),
+                           gko::Array<int>::view(exec,
+                                                 nnz,
+                                                 a.ReadWriteJ(on_device)),
+                           gko::Array<int>::view(exec, a.Height() + 1,
+                                                 a.ReadWriteI(on_device)));
+
+      using ilu_fact_type = gko::factorization::ParIlu<double, int>;
+      std::shared_ptr<ilu_fact_type::Factory> fact_factory = std::move(
+                                                                ilu_fact_type::build()
+                                                                .with_iterations(par_ilu_its)
+                                                                .on(exec));
+
+      if (trisolve_type == "isai")
+      {
+
+         using l_solver_type = gko::preconditioner::LowerIsai<>;
+         using u_solver_type = gko::preconditioner::UpperIsai<>;
+
+         std::shared_ptr<l_solver_type::Factory> l_solver_factory = std::move(
+                                                                       l_solver_type::build()
+                                                                       .with_sparsity_power(sparsity_power)
+                                                                       .on(exec));
+         std::shared_ptr<u_solver_type::Factory> u_solver_factory = std::move(
+                                                                       u_solver_type::build()
+                                                                       .with_sparsity_power(sparsity_power)
+                                                                       .on(exec));
+
+         gko_precond_factory_ = gko::preconditioner::Ilu<l_solver_type,
+         u_solver_type>::build()
+         .with_factorization_factory(fact_factory)
+         .with_l_solver_factory(l_solver_factory)
+         .with_u_solver_factory(u_solver_factory)
+         .on(exec);
+      }
+      else
+      {
+
+         gko_precond_factory_ = gko::preconditioner::Ilu<>::build()
+                                .with_factorization_factory(fact_factory)
+                                .on(exec);
+      }
+
+      gko_precond_ = gko_precond_factory_.get()->generate(
+                        gko::give(gko_sparse));
+   }
+};
+
+class GinkgoCuIcPreconditioner : public GinkgoPreconditionerBase
+{
+public:
+   GinkgoCuIcPreconditioner(std::shared_ptr<const gko::Executor> exec,
+                             SparseMatrix &a, const char *trisolve_type = "exact",
+                             int sparsity_power=1,
+                             bool iter_mode=false)
+      : GinkgoPreconditionerBase(exec, a, iter_mode)
+   {
+
+      bool on_device = false;
+      if (exec->get_master() != exec)
+      {
+         on_device = true;
+      }
+
+      using mtx = gko::matrix::Csr<double, int>;
+      const int nnz =  a.GetMemoryData().Capacity();
+      auto gko_sparse = mtx::create(
+                           exec, gko::dim<2>(a.Height(), a.Width()),
+                           gko::Array<double>::view(exec,
+                                                    nnz,
+                                                    a.ReadWriteData(on_device)),
+                           gko::Array<int>::view(exec,
+                                                 nnz,
+                                                 a.ReadWriteJ(on_device)),
+                           gko::Array<int>::view(exec, a.Height() + 1,
+                                                 a.ReadWriteI(on_device)));
+
+      // TEST
+      gko_sparse->sort_by_column_index();
+
+      using ic_fact_type = gko::factorization::Ic<double, int>;
+      std::shared_ptr<ic_fact_type::Factory> fact_factory = std::move(
+                                                                ic_fact_type::build()
+                                                                .with_both_factors(false)
+                                                                .on(exec));
+
+
+      if (trisolve_type == "isai")
+      {
+
+         using l_solver_type = gko::preconditioner::LowerIsai<>;
+
+         std::shared_ptr<l_solver_type::Factory> l_solver_factory = std::move(
+                                                                       l_solver_type::build()
+                                                                       .with_sparsity_power(sparsity_power)
+                                                                       .on(exec));
+
+         gko_precond_factory_ = gko::preconditioner::Ic<l_solver_type>::build()
+         .with_factorization_factory(fact_factory)
+         .with_l_solver_factory(l_solver_factory)
+         .on(exec);
+
+      }
+      else
+      {
+
+         gko_precond_factory_ = gko::preconditioner::Ic<>::build()
+                                .with_factorization_factory(fact_factory)
+                                .on(exec);
+      }
+
+      gko_precond_ = gko_precond_factory_.get()->generate(
+                        gko::give(gko_sparse));
+   }
+
+   GinkgoCuIcPreconditioner(std::shared_ptr<const gko::Executor> exec,
+                             SparseMatrix &a, Array<int> &inv_permutation_indices,
+                             const char *trisolve_type = "exact",
+                             int sparsity_power=1,
+                             bool iter_mode=false)
+      : GinkgoPreconditionerBase(exec, a, inv_permutation_indices, iter_mode)
+   {
+
+      bool on_device = false;
+      if (exec->get_master() != exec)
+      {
+         on_device = true;
+      }
+
+      using mtx = gko::matrix::Csr<double, int>;
+      const int nnz =  a.GetMemoryData().Capacity();
+      auto gko_sparse = mtx::create(
+                           exec, gko::dim<2>(a.Height(), a.Width()),
+                           gko::Array<double>::view(exec,
+                                                    nnz,
+                                                    a.ReadWriteData(on_device)),
+                           gko::Array<int>::view(exec,
+                                                 nnz,
+                                                 a.ReadWriteJ(on_device)),
+                           gko::Array<int>::view(exec, a.Height() + 1,
+                                                 a.ReadWriteI(on_device)));
+
+      // TEST
+      gko_sparse->sort_by_column_index();
+
+      using ic_fact_type = gko::factorization::Ic<double, int>;
+      std::shared_ptr<ic_fact_type::Factory> fact_factory = std::move(
+                                                                ic_fact_type::build()
+                                                                .with_both_factors(false)
+                                                                .on(exec));
+      if (trisolve_type == "isai")
+      {
+
+         using l_solver_type = gko::preconditioner::LowerIsai<>;
+
+         std::shared_ptr<l_solver_type::Factory> l_solver_factory = std::move(
+                                                                       l_solver_type::build()
+                                                                       .with_sparsity_power(sparsity_power)
+                                                                       .on(exec));
+
+         gko_precond_factory_ = gko::preconditioner::Ic<l_solver_type>::build()
+         .with_factorization_factory(fact_factory)
+         .with_l_solver_factory(l_solver_factory)
+         .on(exec);
+      }
+      else
+      {
+
+         gko_precond_factory_ = gko::preconditioner::Ic<>::build()
+                                .with_factorization_factory(fact_factory)
+                                .on(exec);
+      }
+
+      gko_precond_ = gko_precond_factory_.get()->generate(
+                        gko::give(gko_sparse));
+   }
+};
+
+class GinkgoCuIluPreconditioner : public GinkgoPreconditionerBase
+{
+public:
+   GinkgoCuIluPreconditioner(std::shared_ptr<const gko::Executor> exec,
+                             SparseMatrix &a, const char *trisolve_type = "exact",
+                             int sparsity_power=1,
+                             bool iter_mode=false)
+      : GinkgoPreconditionerBase(exec, a, iter_mode)
+   {
+
+      bool on_device = false;
+      if (exec->get_master() != exec)
+      {
+         on_device = true;
+      }
+
+      using mtx = gko::matrix::Csr<double, int>;
+      const int nnz =  a.GetMemoryData().Capacity();
+      auto gko_sparse = mtx::create(
+                           exec, gko::dim<2>(a.Height(), a.Width()),
+                           gko::Array<double>::view(exec,
+                                                    nnz,
+                                                    a.ReadWriteData(on_device)),
+                           gko::Array<int>::view(exec,
+                                                 nnz,
+                                                 a.ReadWriteJ(on_device)),
+                           gko::Array<int>::view(exec, a.Height() + 1,
+                                                 a.ReadWriteI(on_device)));
+
+      // TEST
+      gko_sparse->sort_by_column_index();
+
+      using ilu_fact_type = gko::factorization::Ilu<double, int>;
+      std::shared_ptr<ilu_fact_type::Factory> fact_factory = std::move(
+                                                                ilu_fact_type::build()
+                                                                .on(exec));
+
+
+      if (trisolve_type == "isai")
+      {
+
+         using l_solver_type = gko::preconditioner::LowerIsai<>;
+         using u_solver_type = gko::preconditioner::UpperIsai<>;
+
+         std::shared_ptr<l_solver_type::Factory> l_solver_factory = std::move(
+                                                                       l_solver_type::build()
+                                                                       .with_sparsity_power(sparsity_power)
+                                                                       .on(exec));
+         std::shared_ptr<u_solver_type::Factory> u_solver_factory = std::move(
+                                                                       u_solver_type::build()
+                                                                       .with_sparsity_power(sparsity_power)
+                                                                       .on(exec));
+
+         gko_precond_factory_ = gko::preconditioner::Ilu<l_solver_type,
+         u_solver_type>::build()
+         .with_factorization_factory(fact_factory)
+         .with_l_solver_factory(l_solver_factory)
+         .with_u_solver_factory(u_solver_factory)
+         .on(exec);
+
+      }
+      else
+      {
+
+         gko_precond_factory_ = gko::preconditioner::Ilu<>::build()
+                                .with_factorization_factory(fact_factory)
+                                .on(exec);
+      }
+
+      gko_precond_ = gko_precond_factory_.get()->generate(
+                        gko::give(gko_sparse));
+   }
+
+   GinkgoCuIluPreconditioner(std::shared_ptr<const gko::Executor> exec,
+                             SparseMatrix &a, Array<int> &inv_permutation_indices,
+                             const char *trisolve_type = "exact",
+                             int sparsity_power=1,
+                             bool iter_mode=false)
+      : GinkgoPreconditionerBase(exec, a, inv_permutation_indices, iter_mode)
+   {
+
+      bool on_device = false;
+      if (exec->get_master() != exec)
+      {
+         on_device = true;
+      }
+
+      using mtx = gko::matrix::Csr<double, int>;
+      const int nnz =  a.GetMemoryData().Capacity();
+      auto gko_sparse = mtx::create(
+                           exec, gko::dim<2>(a.Height(), a.Width()),
+                           gko::Array<double>::view(exec,
+                                                    nnz,
+                                                    a.ReadWriteData(on_device)),
+                           gko::Array<int>::view(exec,
+                                                 nnz,
+                                                 a.ReadWriteJ(on_device)),
+                           gko::Array<int>::view(exec, a.Height() + 1,
+                                                 a.ReadWriteI(on_device)));
+
+      // TEST
+      gko_sparse->sort_by_column_index();
+
+      using ilu_fact_type = gko::factorization::Ilu<double, int>;
+      std::shared_ptr<ilu_fact_type::Factory> fact_factory = std::move(
+                                                                ilu_fact_type::build()
+                                                                .on(exec));
+      if (trisolve_type == "isai")
+      {
+
+         using l_solver_type = gko::preconditioner::LowerIsai<>;
+         using u_solver_type = gko::preconditioner::UpperIsai<>;
+
+         std::shared_ptr<l_solver_type::Factory> l_solver_factory = std::move(
+                                                                       l_solver_type::build()
+                                                                       .with_sparsity_power(sparsity_power)
+                                                                       .on(exec));
+         std::shared_ptr<u_solver_type::Factory> u_solver_factory = std::move(
+                                                                       u_solver_type::build()
+                                                                       .with_sparsity_power(sparsity_power)
+                                                                       .on(exec));
+
+         gko_precond_factory_ = gko::preconditioner::Ilu<l_solver_type,
+         u_solver_type>::build()
+         .with_factorization_factory(fact_factory)
+         .with_l_solver_factory(l_solver_factory)
+         .with_u_solver_factory(u_solver_factory)
+         .on(exec);
+      }
+      else
+      {
+
+         gko_precond_factory_ = gko::preconditioner::Ilu<>::build()
+                                .with_factorization_factory(fact_factory)
+                                .on(exec);
+      }
+
+      gko_precond_ = gko_precond_factory_.get()->generate(
+                        gko::give(gko_sparse));
+   }
+};
+
 } // namespace GinkgoWrappers
 
 }

linalg/solvers.hpp

@@ -28,6 +28,8 @@
 namespace mfem
 {
 
+void MinimumDiscardedFillOrdering(SparseMatrix &C, Array<int> &p);
+
 class BilinearForm;
 
 /// Abstract base class for an iterative solver monitor

mesh/mesh.cpp

@@ -3843,34 +3843,18 @@ Mesh::Mesh(Mesh *orig_mesh, int ref_factor, int ref_type)
 
       rfes.GetBdrElementDofs(el, rdofs);
       MFEM_ASSERT(rdofs.Size() == RG.RefPts.Size(), "");
-      if (Dim == 1)
+      const int *c2h_map = (Dim == 1) ? NULL : rfec.GetDofMap(geom);
+      for (int j = 0; j < RG.RefGeoms.Size()/nvert; j++)
       {
-         // Dim == 1 is a special case because the boundary elements are
-         // zero-dimensional points, and therefore don't have a DofMap
-         for (int j = 0; j < RG.RefGeoms.Size()/nvert; j++)
+         Element *elem = NewElement(geom);
+         elem->SetAttribute(attrib);
+         int *v = elem->GetVertices();
+         for (int k = 0; k < nvert; k++)
          {
-            Element *elem = NewElement(geom);
-            elem->SetAttribute(attrib);
-            int *v = elem->GetVertices();
-            v[0] = rdofs[RG.RefGeoms[nvert*j]];
-            AddBdrElement(elem);
-         }
-      }
-      else
-      {
-         const int *c2h_map = rfec.GetDofMap(geom);
-         for (int j = 0; j < RG.RefGeoms.Size()/nvert; j++)
-         {
-            Element *elem = NewElement(geom);
-            elem->SetAttribute(attrib);
-            int *v = elem->GetVertices();
-            for (int k = 0; k < nvert; k++)
-            {
-               int cid = RG.RefGeoms[k+nvert*j]; // local Cartesian index
-               v[k] = rdofs[c2h_map[cid]];
-            }
-            AddBdrElement(elem);
+            int cid = RG.RefGeoms[k+nvert*j]; // local Cartesian index
+            v[k] = rdofs[c2h_map ? c2h_map[cid] : cid];
          }
+         AddBdrElement(elem);
       }
    }
 
@@ -3911,6 +3895,326 @@ Mesh::Mesh(Mesh *orig_mesh, int ref_factor, int ref_type)
    MFEM_ASSERT(CheckBdrElementOrientation(false) == 0, "");
 }
 
+void Mesh::MakeSimplicial(Mesh &orig_mesh, int *vglobal)
+{
+   MFEM_VERIFY(orig_mesh.CheckElementOrientation(false) == 0,
+               "Mesh::MakeSimplicial requires a properly oriented input mesh");
+
+   int dim = orig_mesh.Dimension();
+   int sdim = orig_mesh.SpaceDimension();
+
+   if (dim == 1)
+   {
+      Mesh copy(orig_mesh);
+      Swap(copy, true);
+      return;
+   }
+
+   int nv = orig_mesh.GetNV();
+   int ne = orig_mesh.GetNE();
+   int nbe = orig_mesh.GetNBE();
+
+   static int num_subdivisions[Geometry::NUM_GEOMETRIES];
+   num_subdivisions[Geometry::POINT] = 1;
+   num_subdivisions[Geometry::SEGMENT] = 1;
+   num_subdivisions[Geometry::TRIANGLE] = 1;
+   num_subdivisions[Geometry::TETRAHEDRON] = 1;
+   num_subdivisions[Geometry::SQUARE] = 2;
+   num_subdivisions[Geometry::PRISM] = 3;
+   num_subdivisions[Geometry::CUBE] = 6;
+   // NOTE: some hexes may be subdivided into only 5 tets, so this is an
+   // estimate only. The actual number of created tets may be less, so the
+   // elements array will need to be shrunk after mesh creation.
+   int new_ne = 0, new_nbe = 0;
+   for (int i=0; i<ne; ++i)
+   {
+      new_ne += num_subdivisions[orig_mesh.GetElementBaseGeometry(i)];
+   }
+   for (int i=0; i<nbe; ++i)
+   {
+      new_nbe += num_subdivisions[orig_mesh.GetBdrElementBaseGeometry(i)];
+   }
+
+   InitMesh(dim, sdim, nv, new_ne, new_nbe);
+
+   // Vertices of the new mesh are same as the original mesh
+   NumOfVertices = nv;
+   for (int i=0; i<nv; ++i)
+   {
+      vertices[i].SetCoords(dim, orig_mesh.vertices[i]());
+   }
+
+   // We need a global vertex numbering to identify which diagonals to split
+   // (quad faces are split using the diagonal originating from the smallest
+   // global vertex number). Use the supplied global numbering, if it is
+   // non-NULL, otherwise use the local numbering.
+   Array<int> vglobal_id;
+   if (vglobal == NULL)
+   {
+      vglobal_id.SetSize(nv);
+      for (int i=0; i<nv; ++i) { vglobal_id[i] = i; }
+      vglobal = vglobal_id.GetData();
+   }
+
+   constexpr int nv_tri = 3, nv_quad = 4, nv_tet = 4, nv_prism = 6, nv_hex = 8;
+   constexpr int quad_ntris = 2, prism_ntets = 3;
+   static const int quad_trimap[2][nv_tri*quad_ntris] =
+   {
+      {
+         0, 0,
+         1, 2,
+         2, 3
+      },{
+         0, 1,
+         1, 2,
+         3, 3
+      }
+   };
+   static const int prism_rot[nv_prism*nv_prism] =
+   {
+      0, 1, 2, 3, 4, 5,
+      1, 2, 0, 4, 5, 3,
+      2, 0, 1, 5, 3, 4,
+      3, 5, 4, 0, 2, 1,
+      4, 3, 5, 1, 0, 2,
+      5, 4, 3, 2, 1, 0
+   };
+   static const int prism_f[nv_quad] = {1, 2, 5, 4};
+   static const int prism_tetmaps[2][nv_prism*prism_ntets] =
+   {
+      {
+         0, 0, 0,
+         1, 1, 4,
+         2, 5, 5,
+         5, 4, 3
+      },{
+         0, 0, 0,
+         1, 4, 4,
+         2, 2, 5,
+         4, 5, 3
+      }
+   };
+   static const int hex_rot[nv_hex*nv_hex] =
+   {
+      0, 1, 2, 3, 4, 5, 6, 7,
+      1, 0, 4, 5, 2, 3, 7, 6,
+      2, 1, 5, 6, 3, 0, 4, 7,
+      3, 0, 1, 2, 7, 4, 5, 6,
+      4, 0, 3, 7, 5, 1, 2, 6,
+      5, 1, 0, 4, 6, 2, 3, 7,
+      6, 2, 1, 5, 7, 3, 0, 4,
+      7, 3, 2, 6, 4, 0, 1, 5
+   };
+   static const int hex_f0[nv_quad] = {1, 2, 6, 5};
+   static const int hex_f1[nv_quad] = {2, 3, 7, 6};
+   static const int hex_f2[nv_quad] = {4, 5, 6, 7};
+   static const int num_rot[8] = {0, 1, 2, 0, 0, 2, 1, 0};
+   static const int hex_tetmap0[nv_tet*5] =
+   {
+      0, 0, 0, 0, 2,
+      1, 2, 2, 5, 7,
+      2, 7, 3, 7, 5,
+      5, 5, 7, 4, 6
+   };
+   static const int hex_tetmap1[nv_tet*6] =
+   {
+      0, 0, 1, 0, 0, 1,
+      5, 1, 6, 7, 7, 7,
+      7, 7, 7, 2, 1, 6,
+      4, 5, 5, 3, 2, 2
+   };
+   static const int hex_tetmap2[nv_tet*6] =
+   {
+      0, 0, 0, 0, 0, 0,
+      4, 3, 7, 1, 3, 6,
+      5, 7, 4, 2, 6, 5,
+      6, 6, 6, 5, 2, 2
+   };
+   static const int hex_tetmap3[nv_tet*6] =
+   {
+      0, 0, 0, 0, 1, 1,
+      2, 3, 7, 5, 5, 6,
+      3, 7, 4, 6, 6, 2,
+      6, 6, 6, 4, 0, 0
+   };
+   static const int *hex_tetmaps[4] =
+   {
+      hex_tetmap0, hex_tetmap1, hex_tetmap2, hex_tetmap3
+   };
+
+   auto find_min = [](const int*a, int n) { return std::min_element(a,a+n)-a; };
+
+   for (int i=0; i<ne; ++i)
+   {
+      const int *v = orig_mesh.elements[i]->GetVertices();
+      const int attrib = orig_mesh.GetAttribute(i);
+      const Geometry::Type orig_geom = orig_mesh.GetElementBaseGeometry(i);
+
+      if (num_subdivisions[orig_geom] == 1)
+      {
+         // (num_subdivisions[orig_geom] == 1) implies that the element does
+         // not need to be further split (it is either a segment, triangle,
+         // or tetrahedron), and so it is left unchanged.
+         Element *e = NewElement(orig_geom);
+         e->SetAttribute(attrib);
+         e->SetVertices(v);
+         AddElement(e);
+      }
+      else if (orig_geom == Geometry::SQUARE)
+      {
+         for (int itri=0; itri<quad_ntris; ++itri)
+         {
+            Element *e = NewElement(Geometry::TRIANGLE);
+            e->SetAttribute(attrib);
+            int *v2 = e->GetVertices();
+            for (int iv=0; iv<nv_tri; ++iv)
+            {
+               v2[iv] = v[quad_trimap[0][itri + iv*quad_ntris]];
+            }
+            AddElement(e);
+         }
+      }
+      else if (orig_geom == Geometry::PRISM)
+      {
+         int vg[nv_prism];
+         for (int iv=0; iv<nv_prism; ++iv) { vg[iv] = vglobal[v[iv]]; }
+         // Rotate the vertices of the prism so that the smallest vertex index
+         // is in the first place
+         int irot = find_min(vg, nv_prism);
+         for (int iv=0; iv<nv_prism; ++iv)
+         {
+            int jv = prism_rot[iv + irot*nv_prism];
+            vg[iv] = v[jv];
+         }
+         // Two cases according to which diagonal splits third quad face
+         int q[nv_quad];
+         for (int iv=0; iv<nv_quad; ++iv) { q[iv] = vglobal[vg[prism_f[iv]]]; }
+         int j = find_min(q, nv_quad);
+         const int *tetmap = (j == 0 || j == 2) ? prism_tetmaps[0] : prism_tetmaps[1];
+         for (int itet=0; itet<prism_ntets; ++itet)
+         {
+            Element *e = NewElement(Geometry::TETRAHEDRON);
+            e->SetAttribute(attrib);
+            int *v2 = e->GetVertices();
+            for (int iv=0; iv<nv_tet; ++iv)
+            {
+               v2[iv] = vg[tetmap[itet + iv*prism_ntets]];
+            }
+            AddElement(e);
+         }
+      }
+      else if (orig_geom == Geometry::CUBE)
+      {
+         int vg[nv_hex];
+         for (int iv=0; iv<nv_hex; ++iv) { vg[iv] = vglobal[v[iv]]; }
+
+         // Rotate the vertices of the hex so that the smallest vertex index is
+         // in the first place
+         int irot = find_min(vg, nv_hex);
+         for (int iv=0; iv<nv_hex; ++iv)
+         {
+            int jv = hex_rot[iv + irot*nv_hex];
+            vg[iv] = v[jv];
+         }
+
+         int q[nv_quad];
+         // Bitmask is three binary digits, each digit is 1 if the diagonal of
+         // the corresponding face goes through the 7th vertex, and 0 if not.
+         int bitmask = 0;
+         int j;
+         // First quad face
+         for (int iv=0; iv<nv_quad; ++iv) { q[iv] = vglobal[vg[hex_f0[iv]]]; }
+         j = find_min(q, nv_quad);
+         if (j == 0 || j == 2) { bitmask += 4; }
+         // Second quad face
+         for (int iv=0; iv<nv_quad; ++iv) { q[iv] = vglobal[vg[hex_f1[iv]]]; }
+         j = find_min(q, nv_quad);
+         if (j == 1 || j == 3) { bitmask += 2; }
+         // Third quad face
+         for (int iv=0; iv<nv_quad; ++iv) { q[iv] = vglobal[vg[hex_f2[iv]]]; }
+         j = find_min(q, nv_quad);
+         if (j == 0 || j == 2) { bitmask += 1; }
+
+         // Apply rotations
+         int nrot = num_rot[bitmask];
+         for (int irot=0; irot<nrot; ++irot)
+         {
+            int vtemp;
+            vtemp = vg[1];
+            vg[1] = vg[4];
+            vg[4] = vg[3];
+            vg[3] = vtemp;
+            vtemp = vg[5];
+            vg[5] = vg[7];
+            vg[7] = vg[2];
+            vg[2] = vtemp;
+         }
+
+         // Sum up nonzero bits in bitmask
+         int ndiags = ((bitmask&4) >> 2) + ((bitmask&2) >> 1) + (bitmask&1);
+         int ntets = (ndiags == 0) ? 5 : 6;
+         const int *tetmap = hex_tetmaps[ndiags];
+         for (int itet=0; itet<ntets; ++itet)
+         {
+            Element *e = NewElement(Geometry::TETRAHEDRON);
+            e->SetAttribute(attrib);
+            int *v2 = e->GetVertices();
+            for (int iv=0; iv<nv_tet; ++iv)
+            {
+               v2[iv] = vg[tetmap[itet + iv*ntets]];
+            }
+            AddElement(e);
+         }
+      }
+   }
+   // In 3D, shrink the element array because some hexes have only 5 tets
+   if (dim == 3) { elements.SetSize(NumOfElements); }
+
+   for (int i=0; i<nbe; ++i)
+   {
+      const int *v = orig_mesh.boundary[i]->GetVertices();
+      const int attrib = orig_mesh.GetBdrAttribute(i);
+      const Geometry::Type orig_geom = orig_mesh.GetBdrElementBaseGeometry(i);
+      if (num_subdivisions[orig_geom] == 1)
+      {
+         Element *be = NewElement(orig_geom);
+         be->SetAttribute(attrib);
+         be->SetVertices(v);
+         AddBdrElement(be);
+      }
+      else if (orig_geom == Geometry::SQUARE)
+      {
+         int vg[nv_quad];
+         for (int iv=0; iv<nv_quad; ++iv) { vg[iv] = vglobal[v[iv]]; }
+         // Split quad according the smallest (global) vertex
+         int iv_min = find_min(vg, nv_quad);
+         int isplit = (iv_min == 0 || iv_min == 2) ? 0 : 1;
+         for (int itri=0; itri<quad_ntris; ++itri)
+         {
+            Element *be = NewElement(Geometry::TRIANGLE);
+            be->SetAttribute(attrib);
+            int *v2 = be->GetVertices();
+            for (int iv=0; iv<nv_tri; ++iv)
+            {
+               v2[iv] = v[quad_trimap[isplit][itri + iv*quad_ntris]];
+            }
+            AddBdrElement(be);
+         }
+      }
+      else
+      {
+         MFEM_ABORT("Unreachable");
+      }
+   }
+
+   FinalizeTopology(false);
+   sequence = orig_mesh.GetSequence();
+   last_operation = orig_mesh.last_operation;
+
+   MFEM_ASSERT(CheckElementOrientation(false) == 0, "");
+   MFEM_ASSERT(CheckBdrElementOrientation(false) == 0, "");
+}
+
 void Mesh::KnotInsert(Array<KnotVector *> &kv)
 {
    if (NURBSext == NULL)
@@ -4174,6 +4478,25 @@ void Mesh::EnsureNodes()
    }
 }
 
+void Mesh::SetVerticesFromNodes()
+{
+   if (!Nodes) { return; }
+   for (int iel=0; iel<GetNE(); ++iel)
+   {
+      Geometry::Type geom = GetElementBaseGeometry(iel);
+      const IntegrationRule *ref_verts = Geometries.GetVertices(geom);
+      DenseMatrix node_coords;
+      ElementTransformation *T = GetElementTransformation(iel);
+      Nodes->GetVectorValues(*T, *ref_verts, node_coords);
+      Element *el = GetElement(iel);
+      for (int iv=0; iv<el->GetNVertices(); ++iv)
+      {
+         int v = el->GetVertices()[iv];
+         vertices[v].SetCoords(&node_coords(0,iv));
+      }
+   }
+}
+
 void Mesh::SetNodalGridFunction(GridFunction *nodes, bool make_owner)
 {
    GetNodes(*nodes);

mesh/mesh.hpp

@@ -701,11 +701,11 @@ public:
                              BasisType::GaussLobatto.
 
        The refinement data which can be accessed with GetRefinementTransforms()
-       is set to reflect the performed refinements.
-
-       @note The constructed Mesh is linear, i.e. it does not have nodes. */
+       is set to reflect the performed refinements. */
    Mesh(Mesh *orig_mesh, int ref_factor, int ref_type);
 
+   void MakeSimplicial(Mesh &orig_mesh, int *vglobal=NULL);
+
    /** This is similar to the mesh constructor with the same arguments, but here
        the current mesh is destroyed and another one created based on the data
        stream again given in MFEM, Netgen, or VTK format. If generate_edges = 0
@@ -1137,6 +1137,9 @@ public:
        with straight edges). */
    void EnsureNodes();
 
+    /** Updates the coordinates of the vertices from the node locations. */
+   void SetVerticesFromNodes();
+
    /** Set the curvature of the mesh nodes using the given polynomial degree,
        'order', and optionally: discontinuous or continuous FE space, 'discont',
        new space dimension, 'space_dim' (if != -1), and 'ordering'. */

mesh/pmesh.cpp

@@ -1142,7 +1142,7 @@ ParMesh::ParMesh(ParMesh *orig_mesh, int ref_factor, int ref_type)
          group_sedge.AddColumnsInRow(gr-1, orig_nq*(RG.RefEdges.Size()/2-
                                                     RG.NumBdrEdges));
          // count refined faces
-         group_squad.AddColumnsInRow(gr-1, orig_nq*(RG.RefGeoms.Size()/nvert));
+         group_squad.AddColumnsInRow(gr-1, orig_nq*RG.RefGeoms.Size()/nvert);
       }
    }
 
@@ -1307,6 +1307,144 @@ ParMesh::ParMesh(ParMesh *orig_mesh, int ref_factor, int ref_type)
    }
 }
 
+void ParMesh::MakeSimplicial(ParMesh &orig_mesh)
+{
+   MyComm = orig_mesh.GetComm();
+   NRanks = orig_mesh.GetNRanks();
+   MyRank = orig_mesh.GetMyRank();
+   glob_elem_offset = -1;
+   glob_offset_sequence = -1;
+   gtopo = orig_mesh.gtopo;
+   have_face_nbr_data = false;
+   pncmesh = NULL;
+   meshgen = orig_mesh.meshgen;
+
+   H1_FECollection fec(1, orig_mesh.Dimension());
+   ParFiniteElementSpace fes(&orig_mesh, &fec);
+
+   Array<int> vglobal(orig_mesh.GetNV());
+   for (int iv=0; iv<orig_mesh.GetNV(); ++iv)
+   {
+      vglobal[iv] = fes.GetGlobalTDofNumber(iv);
+   }
+   Mesh::MakeSimplicial(orig_mesh, vglobal);
+
+   // count the number of entries in each row of group_s{vert,edge,face}
+   group_svert.MakeI(GetNGroups()-1); // exclude the local group 0
+   group_sedge.MakeI(GetNGroups()-1);
+   group_stria.MakeI(GetNGroups()-1);
+   group_squad.MakeI(GetNGroups()-1);
+   for (int gr = 1; gr < GetNGroups(); gr++)
+   {
+      group_svert.AddColumnsInRow(gr-1, orig_mesh.GroupNVertices(gr));
+      group_sedge.AddColumnsInRow(gr-1, orig_mesh.GroupNEdges(gr));
+      // Every quad gives an extra edge
+      const int orig_nq = orig_mesh.GroupNQuadrilaterals(gr);
+      group_sedge.AddColumnsInRow(gr-1, orig_nq);
+      // Every quad is subdivided into two triangles
+      group_stria.AddColumnsInRow(gr-1, 2*orig_nq);
+      // Existing triangles remain unchanged
+      const int orig_nt = orig_mesh.GroupNTriangles(gr);
+      group_stria.AddColumnsInRow(gr-1, orig_nt);
+   }
+   group_svert.MakeJ();
+   svert_lvert.Reserve(group_svert.Size_of_connections());
+
+   group_sedge.MakeJ();
+   shared_edges.Reserve(group_sedge.Size_of_connections());
+   sedge_ledge.SetSize(group_sedge.Size_of_connections());
+
+   group_stria.MakeJ();
+   shared_trias.Reserve(group_stria.Size_of_connections());
+   sface_lface.SetSize(shared_trias.Size());
+
+   group_squad.MakeJ();
+
+   constexpr int ntris = 2, nv_tri = 3, nv_quad = 4;
+
+   Array<int> dofs;
+   for (int gr = 1; gr < GetNGroups(); gr++)
+   {
+      // add shared vertices from original shared vertices
+      const int orig_n_verts = orig_mesh.GroupNVertices(gr);
+      for (int j = 0; j < orig_n_verts; j++)
+      {
+         fes.GetVertexDofs(orig_mesh.GroupVertex(gr, j), dofs);
+         group_svert.AddConnection(gr-1, svert_lvert.Append(dofs[0])-1);
+      }
+
+      // add original shared edges
+      const int orig_n_edges = orig_mesh.GroupNEdges(gr);
+      for (int e = 0; e < orig_n_edges; e++)
+      {
+         int iedge, o;
+         orig_mesh.GroupEdge(gr, e, iedge, o);
+         Element *elem = NewElement(Geometry::SEGMENT);
+         Array<int> edge_verts;
+         orig_mesh.GetEdgeVertices(iedge, edge_verts);
+         elem->SetVertices(edge_verts);
+         group_sedge.AddConnection(gr-1, shared_edges.Append(elem)-1);
+      }
+      // add original shared triangles
+      const int orig_nt = orig_mesh.GroupNTriangles(gr);
+      for (int e = 0; e < orig_nt; e++)
+      {
+         int itri, o;
+         orig_mesh.GroupTriangle(gr, e, itri, o);
+         const int *v = orig_mesh.GetFace(itri)->GetVertices();
+         shared_trias.SetSize(shared_trias.Size()+1);
+         int *v2 = shared_trias.Last().v;
+         for (int iv=0; iv<nv_tri; ++iv) { v2[iv] = v[iv]; }
+         group_stria.AddConnection(gr-1, shared_trias.Size()-1);
+      }
+      // add triangles from split quads and add resulting diagonal edge
+      const int orig_nq = orig_mesh.GroupNQuadrilaterals(gr);
+      if (orig_nq > 0)
+      {
+         static const int trimap[12] =
+         {
+            0, 0, 0, 1,
+            1, 2, 1, 2,
+            2, 3, 3, 3
+         };
+         static const int diagmap[4] = { 0, 2, 1, 3 };
+         for (int f = 0; f < orig_nq; ++f)
+         {
+            int iquad, o;
+            orig_mesh.GroupQuadrilateral(gr, f, iquad, o);
+            const int *v = orig_mesh.GetFace(iquad)->GetVertices();
+            // Split quad according the smallest (global) vertex
+            int vg[nv_quad];
+            for (int iv=0; iv<nv_quad; ++iv) { vg[iv] = vglobal[v[iv]]; }
+            int iv_min = std::min_element(vg, vg+nv_quad) - vg;
+            int isplit = (iv_min == 0 || iv_min == 2) ? 0 : 1;
+            // Add diagonal
+            Element *diag = NewElement(Geometry::SEGMENT);
+            int *v_diag = diag->GetVertices();
+            v_diag[0] = v[diagmap[0 + isplit*2]];
+            v_diag[1] = v[diagmap[1 + isplit*2]];
+            group_sedge.AddConnection(gr-1, shared_edges.Append(diag)-1);
+            // Add two new triangles
+            for (int itri=0; itri<ntris; ++itri)
+            {
+               shared_trias.SetSize(shared_trias.Size()+1);
+               int *v2 = shared_trias.Last().v;
+               for (int iv=0; iv<nv_tri; ++iv)
+               {
+                  v2[iv] = v[trimap[itri + isplit*2 + iv*ntris*2]];
+               }
+               group_stria.AddConnection(gr-1, shared_trias.Size()-1);
+            }
+         }
+      }
+   }
+   group_svert.ShiftUpI();
+   group_sedge.ShiftUpI();
+   group_stria.ShiftUpI();
+
+   FinalizeParTopo();
+}
+
 void ParMesh::Finalize(bool refine, bool fix_orientation)
 {
    const int meshgen_save = meshgen; // Mesh::Finalize() may call SetMeshGen()

mesh/pmesh.hpp

@@ -31,10 +31,10 @@ class ParPumiMesh;
 /// Class for parallel meshes
 class ParMesh : public Mesh
 {
-protected:
+public:
    ParMesh() : MyComm(0), NRanks(0), MyRank(-1),
       have_face_nbr_data(false), pncmesh(NULL) {}
-
+protected:
    MPI_Comm MyComm;
    int NRanks, MyRank;
 
@@ -224,11 +224,11 @@ public:
                              BasisType::GaussLobatto.
 
        The refinement data which can be accessed with GetRefinementTransforms()
-       is set to reflect the performed refinements.
-
-       @note The constructed ParMesh is linear, i.e. it does not have nodes. */
+       is set to reflect the performed refinements. */
    ParMesh(ParMesh *orig_mesh, int ref_factor, int ref_type);
 
+   void MakeSimplicial(ParMesh &orig_mesh);
+
    virtual void Finalize(bool refine = false, bool fix_orientation = false);
 
    virtual void SetAttributes();

miniapps/tools/lor-transfer.cpp

@@ -72,6 +72,7 @@ int main(int argc, char *argv[])
    bool vis = true;
    bool useH1 = false;
    bool use_transfer = false;
+   bool simplex = false;
 
    OptionsParser args(argc, argv);
    args.AddOption(&mesh_file, "-m", "--mesh",
@@ -92,6 +93,8 @@ int main(int argc, char *argv[])
    args.AddOption(&use_transfer, "-t", "--use-pointwise-transfer", "-no-t",
                   "--dont-use-pointwise-transfer",
                   "Use pointwise transfer operators instead of L2 projection.");
+   args.AddOption(&simplex, "-s", "--simplex", "-no-s", "--no-simplex",
+                  "Create a simplicial (tri/tet) LOR mesh.");
    args.Parse();
    if (!args.Good())
    {
@@ -106,7 +109,7 @@ int main(int argc, char *argv[])
 
    // Create the low-order refined mesh
    int basis_lor = BasisType::GaussLobatto; // BasisType::ClosedUniform;
-   Mesh mesh_lor(&mesh, lref, basis_lor);
+   Mesh mesh_lor(&mesh, lref, basis_lor, simplex);
 
    // Create spaces
    FiniteElementCollection *fec, *fec_lor;

tests/unit/mesh/test_mesh.cpp

@@ -145,3 +145,39 @@ TEST_CASE("Gecko integration in MFEM", "[Mesh]")
       }
    }
 }
+
+TEST_CASE("MakeSimplicial", "[Mesh]")
+{
+   auto mesh_fname = GENERATE("../../data/star.mesh",
+                              "../../data/inline-quad.mesh",
+                              "../../data/inline-hex.mesh",
+                              "../../data/inline-wedge.mesh",
+                              "../../data/beam-wedge.mesh");
+
+   Mesh orig_mesh(mesh_fname, 1, 1);
+   Mesh simplex_mesh;
+   simplex_mesh.MakeSimplicial(orig_mesh);
+
+   Geometry::Type orig_geom = orig_mesh.GetElementBaseGeometry(0);
+   int factor;
+   switch (orig_geom)
+   {
+      case Geometry::SQUARE: factor = 2; break;
+      case Geometry::PRISM: factor = 3; break;
+      case Geometry::CUBE: factor = 6; break;
+      default: factor = 1;
+   }
+
+   int dim = orig_mesh.Dimension();
+   Geometry::Type simplex_geom
+      = (dim == 2) ? Geometry::TRIANGLE : Geometry::TETRAHEDRON;
+
+   Array<Geometry::Type> geoms;
+   simplex_mesh.GetGeometries(simplex_mesh.Dimension(), geoms);
+
+   REQUIRE(geoms.Size() == 1);
+   REQUIRE(geoms[0] == simplex_geom);
+   // Note: assuming no hex is subdivided into 5 tets. This can happen depending
+   // on the original mesh, but it doesn't happen for these test cases.
+   REQUIRE(simplex_mesh.GetNE() == orig_mesh.GetNE()*factor);
+}