Updating serial version of example 23

Add IMEX functionality to advection-diffusion example (ex23)

.gitignore

@@ -39,9 +39,6 @@ doc/CodeDocumentation
 *.dSYM
 .DS_Store
 
-# Editor files
-.vscode
-
 # Example and miniapp binaries and outputs
 
 examples/ex[1-9]
@@ -88,7 +85,6 @@ examples/gnuplot_ex20p.inp
 examples/ex22*.mesh
 examples/ex22*.sol
 examples/ex22p_*.*
-examples/exSBP*
 
 examples/sundials/ex9
 examples/sundials/ex1[06]

examples/CMakeLists.txt

@@ -28,6 +28,7 @@ list(APPEND ALL_EXE_SRCS
   ex19.cpp
   ex20.cpp
   ex22.cpp
+  ex23.cpp
   )
 
 if (MFEM_USE_MPI)
@@ -53,6 +54,7 @@ if (MFEM_USE_MPI)
     ex19p.cpp
     ex20p.cpp
     ex22p.cpp
+    ex23p.cpp
     )
 endif()
 

examples/ex23.cpp

@@ -0,0 +1,734 @@
+//                               MFEM Example 23
+//
+// Compile with: make ex23
+//
+// Sample runs:
+//    ex23 -m ../data/periodic-segment.mesh -p 0 -s 2 -dt 0.001 -vs 50
+//    ex23 -m ../data/periodic-segment.mesh -p 0 -s 12 -dt 0.01
+//    ex23 -m ../data/periodic-segment.mesh -p 0 -s 22 -dt 0.01
+//    ex23 -m ../data/periodic-segment.mesh -p 0 -s 32 -dt 0.005 -vs 10
+//    ex23 -m ../data/periodic-square.mesh -p 0 -dt 0.01
+//    ex23 -m ../data/periodic-square.mesh -p 0 -s 32 -dt 0.01
+//    ex23 -m ../data/periodic-hexagon.mesh -p 0 -d 0.001 -s 12 -dt 0.02
+//    ex23 -m ../data/periodic-hexagon.mesh -p 0 -d 0.001 -s 32 -dt 0.009 -vs 10
+//    ex23 -m ../data/periodic-square.mesh -p 1 -dt 0.01 -tf 9
+//    ex23 -m ../data/periodic-hexagon.mesh -p 1 -dt 0.01 -tf 9
+//    ex23 -m ../data/amr-quad.mesh -p 1 -dt 0.01 -tf 9 -vs 2
+//    ex23 -m ../data/disc-nurbs.mesh -p 1 -r 3 -dt 0.01 -tf 9
+//    ex23 -m ../data/disc-nurbs.mesh -p 2 -r 3 -dt 0.01 -tf 9
+//    ex23 -m ../data/disc-nurbs.mesh -p 3 -r 3 -dt 0.01 -tf 9 -d 0.02
+//    ex23 -m ../data/periodic-square.mesh -p 3 -r 3 -dt 0.025 -tf 9
+//    ex23 -m ../data/periodic-cube.mesh -p 0 -o 2 -dt 0.025 -tf 8
+//
+// Description:  This example code solves the time-dependent advection-diffusion
+//               equation
+//               du/dt - div(D grad(u)) + v.grad(u) = 0, where
+//               D is a diffusion coefficient,
+//               v is a given fluid velocity, and
+//               u0(x)=u(0,x) is a given initial condition.
+//
+//               The example demonstrates the use of Discontinuous Galerkin (DG)
+//               bilinear forms in MFEM (face integrators), the use of explicit,
+//               implicit, and implicit-explicit ODE time integrators, the
+//               definition of periodic boundary conditions through periodic
+//               meshes, as well as the use of GLVis for persistent
+//               visualization of a time-evolving solution. The saving of
+//               time-dependent data files for external visualization with
+//               VisIt (visit.llnl.gov) is also illustrated.
+//
+//               This example is a merger of examples 9 and 14.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+// Choice for the problem setup. The fluid velocity, initial condition and
+// boundary condition are chosen based on this parameter.
+int problem;
+
+// Velocity coefficient
+void velocity_function(const Vector &x, Vector &v);
+
+// Initial condition
+double u0_function(const Vector &x);
+
+// Mesh bounding box
+Vector bb_min, bb_max;
+
+
+/** A time-dependent operator for the right-hand side of the ODE for use with
+    explicit ODE solvers. The DG weak form of du/dt = div(D grad(u))-v.grad(u) is
+    M du/dt = - S u + K u + b, where M, S, and K are the mass,
+    stiffness, and advection matrices, and b describes sources and the flow on
+    the boundary.
+    This can be written as a general ODE,
+    du/dt = M^{-1} (-S u + K u + b), and this class is used to compute the RHS
+    and perform the solve for du/dt. */
+class EX_Evolution : public TimeDependentOperator
+{
+private:
+   SparseMatrix &M, &S, &K;
+   const Vector &b;
+
+   DSmoother M_prec;
+   CGSolver M_solver;
+
+   mutable Vector z;
+
+   void initA(double dt);
+
+public:
+   EX_Evolution(SparseMatrix &_M, SparseMatrix &_S, SparseMatrix &_K,
+                const Vector &_b);
+
+   virtual void Mult(const Vector &x, Vector &y) const;
+
+   virtual ~EX_Evolution() {}
+};
+
+/** A time-dependent operator for the right-hand side of the ODE for use with
+    implicit ODE solvers. The DG weak form of du/dt = div(D grad(u))-v.grad(u) is
+    [M + dt (S - K)] du/dt = - S u + K u + b, where M, S, and K are the mass,
+    stiffness, and advection matrices, and b describes sources and the flow on
+    the boundary.
+    This can be written as a general ODE,
+    du/dt = A^{-1} (-S u + K u + b) with A = [M + dt (S - K)], and this class is
+    used to perform the fully implicit solve for du/dt. */
+class IM_Evolution : public TimeDependentOperator
+{
+private:
+   SparseMatrix &M, &S, &K;
+   SparseMatrix *A;
+   const Vector &b;
+
+   DSmoother M_prec;
+   CGSolver M_solver;
+
+   DSmoother *A_prec;
+   GMRESSolver *A_solver;
+   double dt;
+
+   mutable Vector z;
+
+   void initA(double dt);
+
+public:
+   IM_Evolution(SparseMatrix &_M, SparseMatrix &_S, SparseMatrix &_K,
+                const Vector &_b);
+
+   virtual void Mult(const Vector &x, Vector &y) const;
+
+   virtual void ImplicitSolve(const double dt, const Vector &x, Vector &y);
+
+   virtual ~IM_Evolution() { delete A_solver; delete A_prec; delete A; }
+};
+
+/** A time-dependent operator for the right-hand side of the ODE for use with
+    IMEX (Implicit-Explicit) ODE solvers. The DG weak form of
+    du/dt = div(D grad(u))-v.grad(u) is
+    [M + dt S] du/dt = - S u + K u + b, where M, S, and K are the mass,
+    stiffness, and advection matrices, and b describes sources and the flow on
+    the boundary.
+    This can be written as a general ODE,
+    du/dt = A^{-1} (-S u + K u + b) with A = [M + dt (S - K)], and this class is
+    used to perform the implicit or explicit solve for du/dt. */
+class IMEX_Evolution : public TimeDependentOperator
+{
+private:
+   SparseMatrix &M, &S, &K;
+   SparseMatrix *A;
+   const Vector &b;
+
+   DSmoother M_prec;
+   CGSolver M_solver;
+
+   DSmoother *A_prec;
+   CGSolver *A_solver;
+   double dt;
+
+   mutable Vector z;
+
+   void initA(double dt);
+
+public:
+   IMEX_Evolution(SparseMatrix &_M, SparseMatrix &_S, SparseMatrix &_K,
+                  const Vector &_b);
+
+   virtual void ExplicitMult(const Vector &x, Vector &y) const;
+   virtual void Mult(const Vector &x, Vector &y) const;
+
+   virtual void ImplicitSolve(const double dt, const Vector &x, Vector &y);
+
+   virtual ~IMEX_Evolution() { delete A_solver; delete A_prec; delete A; }
+};
+
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   problem = 0;
+   const char *mesh_file = "../data/periodic-hexagon.mesh";
+   int ref_levels = 2;
+   int order = 3;
+   int ode_solver_type = 12;
+   double t_final = 10.0;
+   double d_coef = 0.01;
+   double dt = 0.01;
+   double sigma = -1.0;
+   double kappa = -1.0;
+   bool visualization = true;
+   bool visit = false;
+   bool binary = false;
+   int vis_steps = 5;
+
+   int precision = 8;
+   cout.precision(precision);
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&problem, "-p", "--problem",
+                  "Problem setup to use. See options in velocity_function().");
+   args.AddOption(&ref_levels, "-r", "--refine",
+                  "Number of times to refine the mesh uniformly.");
+   args.AddOption(&order, "-o", "--order",
+                  "Order (degree) of the finite elements.");
+   args.AddOption(&ode_solver_type, "-s", "--ode-solver",
+                  "ODE solver: 1 - Forward Euler, 2 - RK2, 3 - RK3 SSP,"
+                  " 4 - RK4, 5 - Generalized Alpha,\n\t"
+                  "11 - Backward Euler, 12 - SDIRK2, 13 - SDIRK3,\n\t"
+                  "22 - Implicit Midpoint, 23 SDIRK23, 24 - SDIRK34,\n\t"
+                  "31 - IMEX BE/FE, 32 - IMEX RK2.");
+   args.AddOption(&t_final, "-tf", "--t-final",
+                  "Final time; start time is 0.");
+   args.AddOption(&dt, "-dt", "--time-step",
+                  "Time step.");
+   args.AddOption(&d_coef, "-d", "--diff-coef",
+                  "Diffusion coefficient.");
+   args.AddOption(&sigma, "-s", "--sigma",
+                  "One of the two DG penalty parameters, typically +1/-1."
+                  " See the documentation of class DGDiffusionIntegrator.");
+   args.AddOption(&kappa, "-k", "--kappa",
+                  "One of the two DG penalty parameters, should be positive."
+                  " Negative values are replaced with (order+1)^2.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&visit, "-visit", "--visit-datafiles", "-no-visit",
+                  "--no-visit-datafiles",
+                  "Save data files for VisIt (visit.llnl.gov) visualization.");
+   args.AddOption(&binary, "-binary", "--binary-datafiles", "-ascii",
+                  "--ascii-datafiles",
+                  "Use binary (Sidre) or ascii format for VisIt data files.");
+   args.AddOption(&vis_steps, "-vs", "--visualization-steps",
+                  "Visualize every n-th timestep.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   if (kappa < 0)
+   {
+      kappa = (order+1)*(order+1);
+   }
+   args.PrintOptions(cout);
+
+   // 2. Define the ODE solver used for time integration. Several explicit
+   //    Runge-Kutta methods are available.
+   ODESolver *ode_solver = NULL;
+   switch (ode_solver_type)
+   {
+      // Explicit methods
+      case 1:  ode_solver = new ForwardEulerSolver; break;
+      case 2:  ode_solver = new RK2Solver(0.5); break; // midpoint method
+      case 3:  ode_solver = new RK3SSPSolver; break;
+      case 4:  ode_solver = new RK4Solver; break;
+      case 5:  ode_solver = new GeneralizedAlphaSolver(0.5); break;
+      // Implicit L-stable methods
+      case 11: ode_solver = new BackwardEulerSolver; break;
+      case 12: ode_solver = new SDIRK23Solver(2); break;
+      case 13: ode_solver = new SDIRK33Solver; break;
+      // Implicit A-stable methods (not L-stable)
+      case 22: ode_solver = new ImplicitMidpointSolver; break;
+      case 23: ode_solver = new SDIRK23Solver; break;
+      case 24: ode_solver = new SDIRK34Solver; break;
+      // Implicit-Explicit methods
+      case 31: ode_solver = new IMEX_BE_FE; break;
+      case 32: ode_solver = new IMEXRK2; break;
+      default:
+         cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
+         return 3;
+   }
+
+   // 3. Read the serial mesh from the given mesh file on all processors. We can
+   //    handle geometrically periodic meshes in this code.
+   Mesh mesh(mesh_file, 1, 1);
+   int dim = mesh.Dimension();
+
+   // 4. Refine the mesh in serial to increase the resolution. In this example
+   //    we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
+   //    a command-line parameter. If the mesh is of NURBS type, we convert it
+   //    to a (piecewise-polynomial) high-order mesh.
+   for (int lev = 0; lev < ref_levels; lev++)
+   {
+      mesh.UniformRefinement();
+   }
+   if (mesh.NURBSext)
+   {
+      mesh.SetCurvature(max(order, 1));
+   }
+   mesh.GetBoundingBox(bb_min, bb_max, max(order, 1));
+
+   // 5. Define the parallel discontinuous DG finite element space on the
+   //    parallel refined mesh of the given polynomial order.
+   DG_FECollection fec(order, dim);
+   FiniteElementSpace fes(&mesh, &fec);
+
+   cout << "Number of unknowns: " << fes.GetVSize() << endl;
+
+   // 6. Set up and assemble the parallel bilinear and linear forms (and the
+   //    parallel hypre matrices) corresponding to the DG discretization. The
+   //    DGTraceIntegrator involves integrals over mesh interior faces.
+   ConstantCoefficient diff_coef(d_coef);
+   VectorFunctionCoefficient velocity(dim, velocity_function);
+   FunctionCoefficient u0(u0_function);
+
+   BilinearForm m(&fes);
+   m.AddDomainIntegrator(new MassIntegrator);
+
+   BilinearForm s(&fes);
+   s.AddDomainIntegrator(new DiffusionIntegrator(diff_coef));
+   s.AddInteriorFaceIntegrator(new DGDiffusionIntegrator(diff_coef, sigma,
+                                                         kappa));
+   s.AddBdrFaceIntegrator(new DGDiffusionIntegrator(diff_coef, sigma, kappa));
+
+   BilinearForm k(&fes);
+   k.AddDomainIntegrator(new ConvectionIntegrator(velocity, -1.0));
+   k.AddInteriorFaceIntegrator(
+      new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
+   k.AddBdrFaceIntegrator(
+      new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
+
+   LinearForm b(&fes);
+   b.AddBdrFaceIntegrator(
+      new DGDirichletLFIntegrator(u0, diff_coef, sigma, kappa));
+
+   int skip_zeros = 0;
+   m.Assemble(skip_zeros);
+   m.Finalize(skip_zeros);
+   s.Assemble(skip_zeros);
+   s.Finalize(skip_zeros);
+   k.Assemble(skip_zeros);
+   k.Finalize(skip_zeros);
+   b.Assemble();
+
+   // 7. Define the initial conditions, save the corresponding grid function to
+   //    a file and (optionally) save data in the VisIt format and initialize
+   //    GLVis visualization.
+   GridFunction u(&fes);
+   u.ProjectCoefficient(u0);
+
+   {
+      ofstream omesh("ex23.mesh");
+      omesh.precision(precision);
+      mesh.Print(omesh);
+      ofstream osol("ex23-init.gf");
+      osol.precision(precision);
+      u.Save(osol);
+   }
+
+   // Create data collection for solution output: either VisItDataCollection for
+   // ascii data files, or SidreDataCollection for binary data files.
+   DataCollection *dc = NULL;
+   if (visit)
+   {
+      if (binary)
+      {
+#ifdef MFEM_USE_SIDRE
+         dc = new SidreDataCollection("Example23", &mesh);
+#else
+         MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
+#endif
+      }
+      else
+      {
+         dc = new VisItDataCollection("Example23", &mesh);
+         dc->SetPrecision(precision);
+      }
+      dc->RegisterField("solution", &u);
+      dc->SetCycle(0);
+      dc->SetTime(0.0);
+      dc->Save();
+   }
+
+   socketstream sout;
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      sout.open(vishost, visport);
+      if (!sout)
+      {
+         cout << "Unable to connect to GLVis server at "
+              << vishost << ':' << visport << endl;
+         visualization = false;
+         cout << "GLVis visualization disabled.\n";
+      }
+      else
+      {
+         sout.precision(precision);
+         sout << "solution\n" << mesh << u;
+         sout << "pause\n";
+         sout << flush;
+         cout << "GLVis visualization paused."
+              << " Press space (in the GLVis window) to resume it.\n";
+      }
+   }
+
+   // 8. Define the time-dependent evolution operator describing the ODE
+   //    right-hand side, and perform time-integration (looping over the time
+   //    iterations, ti, with a time-step dt).
+
+   TimeDependentOperator *adv = NULL;
+   if (ode_solver_type < 10)
+   {
+      adv = new EX_Evolution(m.SpMat(), s.SpMat(), k.SpMat(), b);
+   }
+   else if (ode_solver_type < 30)
+   {
+      adv = new IM_Evolution(m.SpMat(), s.SpMat(), k.SpMat(), b);
+   }
+   else
+   {
+      adv = new IMEX_Evolution(m.SpMat(), s.SpMat(), k.SpMat(), b);
+   }
+
+   double t = 0.0;
+   adv->SetTime(t);
+   ode_solver->Init(*adv);
+
+   int n_steps = (int)ceil(t_final / dt);
+   double dt_real = t_final / n_steps;
+
+   for (int ti = 0; ti < n_steps; )
+   {
+      ode_solver->Step(u, t, dt_real);
+      ti++;
+
+      if (ti % vis_steps == 0 || ti == n_steps)
+      {
+         cout << "time step: " << ti << ", time: " << t << endl;
+
+         if (visualization)
+         {
+            sout << "solution\n" << mesh << u << flush;
+         }
+
+         if (visit)
+         {
+            dc->SetCycle(ti);
+            dc->SetTime(t);
+            dc->Save();
+         }
+      }
+   }
+
+   // 9. Save the final solution in parallel. This output can be viewed later
+   //    using GLVis: "glvis -np <np> -m ex23-mesh -g ex23-final".
+   {
+      ofstream osol("ex23-final.gf");
+      osol.precision(precision);
+      u.Save(osol);
+   }
+
+   // 10. Free the used memory.
+   delete ode_solver;
+   delete adv;
+   delete dc;
+
+   return 0;
+}
+
+
+// Implementation of class EX_Evolution
+EX_Evolution::EX_Evolution(SparseMatrix &_M, SparseMatrix &_S,
+                           SparseMatrix &_K, const Vector &_b)
+   : TimeDependentOperator(_M.Height()),
+     M(_M), S(_S), K(_K), b(_b), z(_M.Height())
+{
+   M_solver.SetPreconditioner(M_prec);
+   M_solver.SetOperator(M);
+
+   M_solver.iterative_mode = false;
+   M_solver.SetRelTol(1e-9);
+   M_solver.SetAbsTol(0.0);
+   M_solver.SetMaxIter(100);
+   M_solver.SetPrintLevel(0);
+}
+
+void EX_Evolution::Mult(const Vector &x, Vector &y) const
+{
+   // y = M^{-1} (-S x + K x + b)
+   K.Mult(x, z);
+   S.AddMult(x, z, -1.0);
+   z += b;
+   M_solver.Mult(z, y);
+}
+
+// Implementation of class IM_Evolution
+IM_Evolution::IM_Evolution(SparseMatrix &_M, SparseMatrix &_S,
+                           SparseMatrix &_K, const Vector &_b)
+   : TimeDependentOperator(_M.Height()),
+     M(_M), S(_S), K(_K), A(NULL), b(_b),
+     A_prec(NULL), A_solver(NULL), dt(-1.0), z(_M.Height())
+{
+   M_solver.SetPreconditioner(M_prec);
+   M_solver.SetOperator(M);
+
+   M_solver.iterative_mode = false;
+   M_solver.SetRelTol(1e-9);
+   M_solver.SetAbsTol(0.0);
+   M_solver.SetMaxIter(100);
+   M_solver.SetPrintLevel(0);
+}
+
+void IM_Evolution::initA(double _dt)
+{
+   if (fabs(dt - _dt) > 1e-4 * _dt)
+   {
+      delete A_solver;
+      delete A_prec;
+      delete A;
+
+      SparseMatrix * SK = Add(1.0, S, -1.0, K);
+      A = Add(1.0, M, _dt, *SK);
+      delete SK;
+      dt = _dt;
+
+      A_prec = new DSmoother(*A);
+      A_solver = new GMRESSolver;
+      A_solver->SetOperator(*A);
+      A_solver->SetPreconditioner(*A_prec);
+
+      A_solver->iterative_mode = false;
+      A_solver->SetRelTol(1e-9);
+      A_solver->SetAbsTol(0.0);
+      A_solver->SetMaxIter(100);
+      A_solver->SetPrintLevel(0);
+   }
+}
+
+void IM_Evolution::Mult(const Vector &x, Vector &y) const
+{
+   // y = M^{-1} (-S x + K x + b)
+   K.Mult(x, z);
+   S.AddMult(x, z, -1.0);
+   z += b;
+   M_solver.Mult(z, y);
+}
+
+void IM_Evolution::ImplicitSolve(const double _dt, const Vector &x, Vector &y)
+{
+   this->initA(_dt);
+
+   // y = (M + dt S - dt K)^{-1} (-S x + K x + b)
+   K.Mult(x, z);
+   S.AddMult(x, z, -1.0);
+   z += b;
+   A_solver->Mult(z, y);
+}
+
+// Implementation of class IMEX_Evolution
+IMEX_Evolution::IMEX_Evolution(SparseMatrix &_M, SparseMatrix &_S,
+                               SparseMatrix &_K, const Vector &_b)
+   : TimeDependentOperator(_M.Height()),
+     M(_M), S(_S), K(_K), A(NULL), b(_b),
+     A_prec(NULL), A_solver(NULL), dt(-1.0), z(_M.Height())
+{
+   M_solver.SetPreconditioner(M_prec);
+   M_solver.SetOperator(M);
+
+   M_solver.iterative_mode = false;
+   M_solver.SetRelTol(1e-9);
+   M_solver.SetAbsTol(0.0);
+   M_solver.SetMaxIter(100);
+   M_solver.SetPrintLevel(0);
+}
+
+void IMEX_Evolution::initA(double _dt)
+{
+   if (fabs(dt - _dt) > 1e-4 * _dt)
+   {
+      delete A_solver;
+      delete A_prec;
+      delete A;
+
+      A = Add(_dt, S, 1.0, M); // A = M + dt * S
+      dt = _dt;
+
+      A_prec = new DSmoother(*A);
+      A_solver = new CGSolver;
+      A_solver->SetOperator(*A);
+      A_solver->SetPreconditioner(*A_prec);
+
+      A_solver->iterative_mode = false;
+      A_solver->SetRelTol(1e-9);
+      A_solver->SetAbsTol(0.0);
+      A_solver->SetMaxIter(100);
+      A_solver->SetPrintLevel(0);
+   }
+}
+
+void IMEX_Evolution::Mult(const Vector &x, Vector &y) const
+{
+   // y = M^{-1} (-S x + K x + b)
+   K.Mult(x, z);
+   S.AddMult(x, z, -1.0);
+   z += b;
+   M_solver.Mult(z, y);
+}
+
+void IMEX_Evolution::ExplicitMult(const Vector &x, Vector &y) const
+{
+   // y = M^{-1} (K x + b)
+   K.Mult(x, z);
+   z += b;
+   M_solver.Mult(z, y);
+}
+
+void IMEX_Evolution::ImplicitSolve(const double _dt, const Vector &x, Vector &y)
+{
+   this->initA(_dt);
+   // y = (M + dt S)^{-1} (-S x + b)
+   S.Mult(x, z);
+   z *= -1.0;
+   z += b;
+   A_solver->Mult(z, y);
+}
+
+// Velocity coefficient
+void velocity_function(const Vector &x, Vector &v)
+{
+   int dim = x.Size();
+
+   // map to the reference [-1,1] domain
+   Vector X(dim);
+   for (int i = 0; i < dim; i++)
+   {
+      double center = (bb_min[i] + bb_max[i]) * 0.5;
+      X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
+   }
+
+   switch (problem)
+   {
+      case 0:
+      {
+         // Translations in 1D, 2D, and 3D
+         switch (dim)
+         {
+            case 1: v(0) = 1.0; break;
+            case 2: v(0) = sqrt(2./3.); v(1) = sqrt(1./3.); break;
+            case 3: v(0) = sqrt(3./6.); v(1) = sqrt(2./6.); v(2) = sqrt(1./6.);
+               break;
+         }
+         break;
+      }
+      case 1:
+      case 2:
+      {
+         // Clockwise rotation in 2D around the origin
+         const double w = M_PI/2;
+         switch (dim)
+         {
+            case 1: v(0) = 1.0; break;
+            case 2: v(0) = w*X(1); v(1) = -w*X(0); break;
+            case 3: v(0) = w*X(1); v(1) = -w*X(0); v(2) = 0.0; break;
+         }
+         break;
+      }
+      case 3:
+      {
+         // Clockwise twisting rotation in 2D around the origin
+         const double w = M_PI/2;
+         double d = max((X(0)+1.)*(1.-X(0)),0.) * max((X(1)+1.)*(1.-X(1)),0.);
+         d = d*d;
+         switch (dim)
+         {
+            case 1: v(0) = 1.0; break;
+            case 2: v(0) = d*w*X(1); v(1) = -d*w*X(0); break;
+            case 3: v(0) = d*w*X(1); v(1) = -d*w*X(0); v(2) = 0.0; break;
+         }
+         break;
+      }
+   }
+}
+
+// Initial condition
+double u0_function(const Vector &x)
+{
+   int dim = x.Size();
+
+   // map to the reference [-1,1] domain
+   Vector X(dim);
+   for (int i = 0; i < dim; i++)
+   {
+      double center = (bb_min[i] + bb_max[i]) * 0.5;
+      X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
+   }
+
+   switch (problem)
+   {
+      case 0:
+      case 1:
+      {
+         switch (dim)
+         {
+            case 1:
+               return exp(-40.*pow(X(0)-0.5,2));
+            case 2:
+            case 3:
+            {
+               double rx = 0.45, ry = 0.25, cx = 0., cy = -0.2, w = 10.;
+               if (dim == 3)
+               {
+                  const double s = (1. + 0.25*cos(2*M_PI*X(2)));
+                  rx *= s;
+                  ry *= s;
+               }
+               return ( erfc(w*(X(0)-cx-rx))*erfc(-w*(X(0)-cx+rx)) *
+                        erfc(w*(X(1)-cy-ry))*erfc(-w*(X(1)-cy+ry)) )/16;
+            }
+         }
+      }
+      case 2:
+      {
+         double x_ = X(0), y_ = X(1), rho, phi;
+         rho = hypot(x_, y_);
+         phi = atan2(y_, x_);
+         return pow(sin(M_PI*rho),2)*sin(3*phi);
+      }
+      case 3:
+      {
+         const double f = M_PI;
+         return sin(f*X(0))*sin(f*X(1));
+      }
+   }
+   return 0.0;
+}
+
+// Inflow boundary condition (zero for the problems considered in this example)
+double inflow_function(const Vector &x)
+{
+   switch (problem)
+   {
+      case 0:
+      case 1:
+      case 2:
+      case 3: return 0.0;
+   }
+   return 0.0;
+}

examples/ex23p.cpp

@@ -0,0 +1,797 @@
+//                       MFEM Example 23 - Parallel Version
+//
+// Compile with: make ex23p
+//
+// Sample runs:
+//    mpirun -np 4 ex23p -m ../data/periodic-segment.mesh -p 0 -s 2 -dt 0.001 -vs 50
+//    mpirun -np 4 ex23p -m ../data/periodic-segment.mesh -p 0 -s 12 -dt 0.01
+//    mpirun -np 4 ex23p -m ../data/periodic-segment.mesh -p 0 -s 22 -dt 0.01
+//    mpirun -np 4 ex23p -m ../data/periodic-segment.mesh -p 0 -s 32 -dt 0.005 -vs 10
+//    mpirun -np 4 ex23p -m ../data/periodic-square.mesh -p 0 -dt 0.01
+//    mpirun -np 4 ex23p -m ../data/periodic-square.mesh -p 0 -s 32 -dt 0.01
+//    mpirun -np 4 ex23p -m ../data/periodic-hexagon.mesh -p 0 -d 0.001 -s 12 -dt 0.02
+//    mpirun -np 4 ex23p -m ../data/periodic-hexagon.mesh -p 0 -d 0.001 -s 32 -dt 0.009 -vs 10
+//    mpirun -np 4 ex23p -m ../data/periodic-square.mesh -p 1 -dt 0.01 -tf 9
+//    mpirun -np 4 ex23p -m ../data/periodic-hexagon.mesh -p 1 -dt 0.01 -tf 9
+//    mpirun -np 4 ex23p -m ../data/amr-quad.mesh -p 1 -dt 0.01 -tf 9 -vs 2
+//    mpirun -np 4 ex23p -m ../data/disc-nurbs.mesh -p 1 -rp 1 -dt 0.01 -tf 9
+//    mpirun -np 4 ex23p -m ../data/disc-nurbs.mesh -p 2 -rp 1 -dt 0.01 -tf 9
+//    mpirun -np 4 ex23p -m ../data/disc-nurbs.mesh -p 3 -rp 1 -dt 0.01 -tf 9 -d 0.02
+//    mpirun -np 4 ex23p -m ../data/periodic-square.mesh -p 3 -rp 1 -dt 0.025 -tf 9
+//    mpirun -np 4 ex23p -m ../data/periodic-cube.mesh -p 0 -o 2 -dt 0.025 -tf 8
+//
+// Description:  This example code solves the time-dependent advection-diffusion
+//               equation
+//               du/dt - div(D grad(u)) + v.grad(u) = 0, where
+//               D is a diffusion coefficient,
+//               v is a given fluid velocity, and
+//               u0(x)=u(0,x) is a given initial condition.
+//
+//               The example demonstrates the use of Discontinuous Galerkin (DG)
+//               bilinear forms in MFEM (face integrators), the use of explicit,
+//               implicit, and implicit-explicit ODE time integrators, the
+//               definition of periodic boundary conditions through periodic
+//               meshes, as well as the use of GLVis for persistent
+//               visualization of a time-evolving solution. The saving of
+//               time-dependent data files for external visualization with
+//               VisIt (visit.llnl.gov) is also illustrated.
+//
+//               This example is a merger of examples 9 and 14.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+// Choice for the problem setup. The fluid velocity, initial condition and
+// boundary condition are chosen based on this parameter.
+int problem;
+
+// Velocity coefficient
+void velocity_function(const Vector &x, Vector &v);
+
+// Initial condition
+double u0_function(const Vector &x);
+
+// Mesh bounding box
+Vector bb_min, bb_max;
+
+/** A time-dependent operator for the right-hand side of the ODE for use with
+    explicit ODE solvers. The DG weak form of du/dt = div(D grad(u))-v.grad(u) is
+    M du/dt = - S u + K u + b, where M, S, and K are the mass,
+    stiffness, and advection matrices, and b describes sources and the flow on
+    the boundary.
+    This can be written as a general ODE,
+    du/dt = M^{-1} (-S u + K u + b), and this class is used to compute the RHS
+    and perform the solve for du/dt. */
+class EX_Evolution : public TimeDependentOperator
+{
+private:
+   HypreParMatrix &M, &S, &K;
+   const Vector &b;
+
+   HypreSmoother M_prec;
+   CGSolver M_solver;
+
+   mutable Vector z;
+
+   void initA(double dt);
+
+public:
+   EX_Evolution(HypreParMatrix &_M, HypreParMatrix &_S, HypreParMatrix &_K,
+                const Vector &_b);
+
+   virtual void Mult(const Vector &x, Vector &y) const;
+
+   virtual ~EX_Evolution() {}
+};
+
+/** A time-dependent operator for the right-hand side of the ODE for use with
+    implicit ODE solvers. The DG weak form of du/dt = div(D grad(u))-v.grad(u) is
+    [M + dt (S - K)] du/dt = - S u + K u + b, where M, S, and K are the mass,
+    stiffness, and advection matrices, and b describes sources and the flow on
+    the boundary.
+    This can be written as a general ODE,
+    du/dt = A^{-1} (-S u + K u + b) with A = [M + dt (S - K)], and this class is
+    used to perform the fully implicit solve for du/dt. */
+class IM_Evolution : public TimeDependentOperator
+{
+private:
+   HypreParMatrix &M, &S, &K;
+   HypreParMatrix *A;
+   const Vector &b;
+
+   HypreSmoother M_prec;
+   CGSolver M_solver;
+
+   HypreBoomerAMG *A_prec;
+   GMRESSolver *A_solver;
+   double dt;
+
+   mutable Vector z;
+
+   void initA(double dt);
+
+public:
+   IM_Evolution(HypreParMatrix &_M, HypreParMatrix &_S, HypreParMatrix &_K,
+                const Vector &_b);
+
+   virtual void Mult(const Vector &x, Vector &y) const;
+
+   virtual void ImplicitSolve(const double dt, const Vector &x, Vector &y);
+
+   virtual ~IM_Evolution() { delete A_solver; delete A_prec; delete A; }
+};
+
+/** A time-dependent operator for the right-hand side of the ODE for use with
+    IMEX (Implicit-Explicit) ODE solvers. The DG weak form of
+    du/dt = div(D grad(u))-v.grad(u) is
+    [M + dt S] du/dt = - S u + K u + b, where M, S, and K are the mass,
+    stiffness, and advection matrices, and b describes sources and the flow on
+    the boundary.
+    This can be written as a general ODE,
+    du/dt = A^{-1} (-S u + K u + b) with A = [M + dt (S - K)], and this class is
+    used to perform the implicit or explicit solve for du/dt. */
+class IMEX_Evolution : public TimeDependentOperator
+{
+private:
+   HypreParMatrix &M, &S, &K;
+   HypreParMatrix *A;
+   const Vector &b;
+
+   HypreSmoother M_prec;
+   CGSolver M_solver;
+
+   HypreBoomerAMG *A_prec;
+   CGSolver *A_solver;
+   double dt;
+
+   mutable Vector z;
+
+   void initA(double dt);
+
+public:
+   IMEX_Evolution(HypreParMatrix &_M, HypreParMatrix &_S, HypreParMatrix &_K,
+                  const Vector &_b);
+
+   virtual void ExplicitMult(const Vector &x, Vector &y) const;
+   virtual void Mult(const Vector &x, Vector &y) const;
+
+   virtual void ImplicitSolve(const double dt, const Vector &x, Vector &y);
+
+   virtual ~IMEX_Evolution() { delete A_solver; delete A_prec; delete A; }
+};
+
+
+int main(int argc, char *argv[])
+{
+   // 1. Initialize MPI.
+   int num_procs, myid;
+   MPI_Init(&argc, &argv);
+   MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
+   MPI_Comm_rank(MPI_COMM_WORLD, &myid);
+
+   // 2. Parse command-line options.
+   problem = 0;
+   const char *mesh_file = "../data/periodic-hexagon.mesh";
+   int ser_ref_levels = 2;
+   int par_ref_levels = 0;
+   int order = 3;
+   int ode_solver_type = 12;
+   double t_final = 10.0;
+   double d_coef = 0.01;
+   double dt = 0.01;
+   double sigma = -1.0;
+   double kappa = -1.0;
+   bool visualization = true;
+   bool visit = false;
+   bool binary = false;
+   int vis_steps = 5;
+
+   int precision = 8;
+   cout.precision(precision);
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&problem, "-p", "--problem",
+                  "Problem setup to use. See options in velocity_function().");
+   args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
+                  "Number of times to refine the mesh uniformly in serial.");
+   args.AddOption(&par_ref_levels, "-rp", "--refine-parallel",
+                  "Number of times to refine the mesh uniformly in parallel.");
+   args.AddOption(&order, "-o", "--order",
+                  "Order (degree) of the finite elements.");
+   args.AddOption(&ode_solver_type, "-s", "--ode-solver",
+                  "ODE solver: 1 - Forward Euler, 2 - RK2, 3 - RK3 SSP,"
+                  " 4 - RK4, 5 - Generalized Alpha,\n\t"
+                  "11 - Backward Euler, 12 - SDIRK2, 13 - SDIRK3,\n\t"
+                  "22 - Implicit Midpoint, 23 SDIRK23, 24 - SDIRK34,\n\t"
+                  "31 - IMEX BE/FE, 32 - IMEX RK2.");
+   args.AddOption(&t_final, "-tf", "--t-final",
+                  "Final time; start time is 0.");
+   args.AddOption(&dt, "-dt", "--time-step",
+                  "Time step.");
+   args.AddOption(&d_coef, "-d", "--diff-coef",
+                  "Diffusion coefficient.");
+   args.AddOption(&sigma, "-s", "--sigma",
+                  "One of the two DG penalty parameters, typically +1/-1."
+                  " See the documentation of class DGDiffusionIntegrator.");
+   args.AddOption(&kappa, "-k", "--kappa",
+                  "One of the two DG penalty parameters, should be positive."
+                  " Negative values are replaced with (order+1)^2.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&visit, "-visit", "--visit-datafiles", "-no-visit",
+                  "--no-visit-datafiles",
+                  "Save data files for VisIt (visit.llnl.gov) visualization.");
+   args.AddOption(&binary, "-binary", "--binary-datafiles", "-ascii",
+                  "--ascii-datafiles",
+                  "Use binary (Sidre) or ascii format for VisIt data files.");
+   args.AddOption(&vis_steps, "-vs", "--visualization-steps",
+                  "Visualize every n-th timestep.");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      MPI_Finalize();
+      return 1;
+   }
+   if (kappa < 0)
+   {
+      kappa = (order+1)*(order+1);
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+
+   // 3. Define the ODE solver used for time integration. Several explicit,
+   //    implicitit, and implicit-explicit Runge-Kutta methods are available.
+   ODESolver *ode_solver = NULL;
+
+   switch (ode_solver_type)
+   {
+      // Explicit methods
+      case 1:  ode_solver = new ForwardEulerSolver; break;
+      case 2:  ode_solver = new RK2Solver(0.5); break; // midpoint method
+      case 3:  ode_solver = new RK3SSPSolver; break;
+      case 4:  ode_solver = new RK4Solver; break;
+      case 5:  ode_solver = new GeneralizedAlphaSolver(0.5); break;
+      // Implicit L-stable methods
+      case 11: ode_solver = new BackwardEulerSolver; break;
+      case 12: ode_solver = new SDIRK23Solver(2); break;
+      case 13: ode_solver = new SDIRK33Solver; break;
+      // Implicit A-stable methods (not L-stable)
+      case 22: ode_solver = new ImplicitMidpointSolver; break;
+      case 23: ode_solver = new SDIRK23Solver; break;
+      case 24: ode_solver = new SDIRK34Solver; break;
+      // Implicit-Explicit methods
+      case 31: ode_solver = new IMEX_BE_FE; break;
+      case 32: ode_solver = new IMEXRK2; break;
+      default:
+         cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
+         return 3;
+   }
+
+   // 4. Read the serial mesh from the given mesh file on all processors. We can
+   //    handle geometrically periodic meshes in this code.
+   Mesh *mesh = new Mesh(mesh_file, 1, 1);
+   int dim = mesh->Dimension();
+
+   // 5. Refine the mesh in serial to increase the resolution. In this example
+   //    we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
+   //    a command-line parameter. If the mesh is of NURBS type, we convert it
+   //    to a (piecewise-polynomial) high-order mesh.
+   for (int lev = 0; lev < ser_ref_levels; lev++)
+   {
+      mesh->UniformRefinement();
+   }
+   if (mesh->NURBSext)
+   {
+      mesh->SetCurvature(max(order, 1));
+   }
+   mesh->GetBoundingBox(bb_min, bb_max, max(order, 1));
+
+   // 6. Define the parallel mesh by a partitioning of the serial mesh. Refine
+   //    this mesh further in parallel to increase the resolution. Once the
+   //    parallel mesh is defined, the serial mesh can be deleted.
+   ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
+   delete mesh;
+   for (int lev = 0; lev < par_ref_levels; lev++)
+   {
+      pmesh->UniformRefinement();
+   }
+
+   // 7. Define the parallel discontinuous DG finite element space on the
+   //    parallel refined mesh of the given polynomial order.
+   DG_FECollection fec(order, dim);
+   ParFiniteElementSpace *fes = new ParFiniteElementSpace(pmesh, &fec);
+
+   HYPRE_Int global_vSize = fes->GlobalTrueVSize();
+   if (myid == 0)
+   {
+      cout << "Number of unknowns: " << global_vSize << endl;
+   }
+
+   // 8. Set up and assemble the parallel bilinear and linear forms (and the
+   //    parallel hypre matrices) corresponding to the DG discretization. The
+   //    DGTraceIntegrator involves integrals over mesh interior faces.
+   ConstantCoefficient diff_coef(d_coef);
+   VectorFunctionCoefficient velocity(dim, velocity_function);
+   FunctionCoefficient u0(u0_function);
+
+   ParBilinearForm *m = new ParBilinearForm(fes);
+   m->AddDomainIntegrator(new MassIntegrator);
+
+   ParBilinearForm *s = new ParBilinearForm(fes);
+   s->AddDomainIntegrator(new DiffusionIntegrator(diff_coef));
+   s->AddInteriorFaceIntegrator(new DGDiffusionIntegrator(diff_coef, sigma,
+                                                          kappa));
+   s->AddBdrFaceIntegrator(new DGDiffusionIntegrator(diff_coef, sigma, kappa));
+
+   ParBilinearForm *k = new ParBilinearForm(fes);
+   k->AddDomainIntegrator(new ConvectionIntegrator(velocity, -1.0));
+   k->AddInteriorFaceIntegrator(
+      new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
+   k->AddBdrFaceIntegrator(
+      new TransposeIntegrator(new DGTraceIntegrator(velocity, 1.0, -0.5)));
+
+   ParLinearForm *b = new ParLinearForm(fes);
+   b->AddBdrFaceIntegrator(
+      new DGDirichletLFIntegrator(u0, diff_coef, sigma, kappa));
+
+   int skip_zeros = 0;
+   m->Assemble(skip_zeros);
+   m->Finalize(skip_zeros);
+   s->Assemble(skip_zeros);
+   s->Finalize(skip_zeros);
+   k->Assemble(skip_zeros);
+   k->Finalize(skip_zeros);
+   b->Assemble();
+
+   HypreParMatrix *M = m->ParallelAssemble();
+   HypreParMatrix *S = s->ParallelAssemble();
+   HypreParMatrix *K = k->ParallelAssemble();
+   HypreParVector *B = b->ParallelAssemble();
+
+   // 9. Define the initial conditions, save the corresponding grid function to
+   //    a file and (optionally) save data in the VisIt format and initialize
+   //    GLVis visualization.
+   ParGridFunction *u = new ParGridFunction(fes);
+   u->ProjectCoefficient(u0);
+   HypreParVector *U = u->GetTrueDofs();
+
+   {
+      ostringstream mesh_name, sol_name;
+      mesh_name << "ex23-mesh." << setfill('0') << setw(6) << myid;
+      sol_name << "ex23-init." << setfill('0') << setw(6) << myid;
+      ofstream omesh(mesh_name.str().c_str());
+      omesh.precision(precision);
+      pmesh->Print(omesh);
+      ofstream osol(sol_name.str().c_str());
+      osol.precision(precision);
+      u->Save(osol);
+   }
+
+   // Create data collection for solution output: either VisItDataCollection for
+   // ascii data files, or SidreDataCollection for binary data files.
+   DataCollection *dc = NULL;
+   if (visit)
+   {
+      if (binary)
+      {
+#ifdef MFEM_USE_SIDRE
+         dc = new SidreDataCollection("Example23-Parallel", pmesh);
+#else
+         MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
+#endif
+      }
+      else
+      {
+         dc = new VisItDataCollection("Example23-Parallel", pmesh);
+         dc->SetPrecision(precision);
+         // To save the mesh using MFEM's parallel mesh format:
+         // dc->SetFormat(DataCollection::PARALLEL_FORMAT);
+      }
+      dc->RegisterField("solution", u);
+      dc->SetCycle(0);
+      dc->SetTime(0.0);
+      dc->Save();
+   }
+
+   socketstream sout;
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      sout.open(vishost, visport);
+      if (!sout)
+      {
+         if (myid == 0)
+            cout << "Unable to connect to GLVis server at "
+                 << vishost << ':' << visport << endl;
+         visualization = false;
+         if (myid == 0)
+         {
+            cout << "GLVis visualization disabled.\n";
+         }
+      }
+      else
+      {
+         sout << "parallel " << num_procs << " " << myid << "\n";
+         sout.precision(precision);
+         sout << "solution\n" << *pmesh << *u;
+         sout << "pause\n";
+         sout << flush;
+         if (myid == 0)
+            cout << "GLVis visualization paused."
+                 << " Press space (in the GLVis window) to resume it.\n";
+      }
+   }
+
+   // 10. Define the time-dependent evolution operator describing the ODE
+   //     right-hand side, and perform time-integration (looping over the time
+   //     iterations, ti, with a time-step dt).
+
+   TimeDependentOperator *adv = NULL;
+   if (ode_solver_type < 10)
+   {
+      adv = new EX_Evolution(*M, *S, *K, *B);
+   }
+   else if (ode_solver_type < 30)
+   {
+      adv = new IM_Evolution(*M, *S, *K, *B);
+   }
+   else
+   {
+      adv = new IMEX_Evolution(*M, *S, *K, *B);
+   }
+
+   double t = 0.0;
+   adv->SetTime(t);
+   ode_solver->Init(*adv);
+
+   int n_steps = (int)ceil(t_final / dt);
+   double dt_real = t_final / n_steps;
+
+   for (int ti = 0; ti < n_steps; )
+   {
+      ode_solver->Step(*U, t, dt_real);
+      ti++;
+
+      if (ti % vis_steps == 0 || ti == n_steps)
+      {
+         if (myid == 0)
+         {
+            cout << "time step: " << ti << ", time: " << t << endl;
+         }
+
+         // 11. Extract the parallel grid function corresponding to the finite
+         //     element approximation U (the local solution on each processor).
+         *u = *U;
+
+         if (visualization)
+         {
+            sout << "parallel " << num_procs << " " << myid << "\n";
+            sout << "solution\n" << *pmesh << *u << flush;
+         }
+
+         if (visit)
+         {
+            dc->SetCycle(ti);
+            dc->SetTime(t);
+            dc->Save();
+         }
+      }
+   }
+
+   // 12. Save the final solution in parallel. This output can be viewed later
+   //     using GLVis: "glvis -np <np> -m ex23-mesh -g ex23-final".
+   {
+      *u = *U;
+      ostringstream sol_name;
+      sol_name << "ex23-final." << setfill('0') << setw(6) << myid;
+      ofstream osol(sol_name.str().c_str());
+      osol.precision(precision);
+      u->Save(osol);
+   }
+
+   // 13. Free the used memory.
+   delete U;
+   delete u;
+   delete B;
+   delete b;
+   delete K;
+   delete k;
+   delete S;
+   delete s;
+   delete M;
+   delete m;
+   delete fes;
+   delete pmesh;
+   delete ode_solver;
+   delete adv;
+   delete dc;
+
+   MPI_Finalize();
+   return 0;
+}
+
+
+// Implementation of class EX_Evolution
+EX_Evolution::EX_Evolution(HypreParMatrix &_M, HypreParMatrix &_S,
+                           HypreParMatrix &_K, const Vector &_b)
+   : TimeDependentOperator(_M.Height()),
+     M(_M), S(_S), K(_K), b(_b),
+     M_prec(M), M_solver(M.GetComm()), z(M.Height())
+{
+   M_prec.SetType(HypreSmoother::Jacobi);
+   M_solver.SetPreconditioner(M_prec);
+   M_solver.SetOperator(M);
+
+   M_solver.iterative_mode = false;
+   M_solver.SetRelTol(1e-9);
+   M_solver.SetAbsTol(0.0);
+   M_solver.SetMaxIter(100);
+   M_solver.SetPrintLevel(0);
+}
+
+void EX_Evolution::Mult(const Vector &x, Vector &y) const
+{
+   // y = M^{-1} (-S x + K x + b)
+   S.Mult(-1.0, x, 0.0, z);
+   K.Mult(1.0, x, 1.0, z);
+   z += b;
+   M_solver.Mult(z, y);
+}
+
+// Implementation of class IM_Evolution
+IM_Evolution::IM_Evolution(HypreParMatrix &_M, HypreParMatrix &_S,
+                           HypreParMatrix &_K, const Vector &_b)
+   : TimeDependentOperator(_M.Height()),
+     M(_M), S(_S), K(_K), A(NULL), b(_b),
+     M_prec(M), M_solver(M.GetComm()),
+     A_prec(NULL), A_solver(NULL), dt(-1.0), z(M.Height())
+{
+   M_prec.SetType(HypreSmoother::Jacobi);
+   M_solver.SetPreconditioner(M_prec);
+   M_solver.SetOperator(M);
+
+   M_solver.iterative_mode = false;
+   M_solver.SetRelTol(1e-9);
+   M_solver.SetAbsTol(0.0);
+   M_solver.SetMaxIter(100);
+   M_solver.SetPrintLevel(0);
+}
+
+void IM_Evolution::initA(double _dt)
+{
+   if (fabs(dt - _dt) > 1e-4 * _dt)
+   {
+      delete A_solver;
+      delete A_prec;
+      delete A;
+
+      HypreParMatrix * SK = Add(1.0, S, -1.0, K); // SK = S - K
+      A = Add(_dt, *SK, 1.0, M);                  // A = M + dt * (S - K)
+      delete SK;
+      dt = _dt;
+
+      A_prec = new HypreBoomerAMG(*A);
+      A_solver = new GMRESSolver(A->GetComm());
+      A_solver->SetOperator(*A);
+      A_solver->SetPreconditioner(*A_prec);
+
+      A_solver->iterative_mode = false;
+      A_solver->SetRelTol(1e-9);
+      A_solver->SetAbsTol(0.0);
+      A_solver->SetMaxIter(100);
+      A_solver->SetPrintLevel(0);
+   }
+}
+
+void IM_Evolution::Mult(const Vector &x, Vector &y) const
+{
+   // y = M^{-1} (-S x + K x + b)
+   S.Mult(-1.0, x, 0.0, z);
+   K.Mult(1.0, x, 1.0, z);
+   z += b;
+   M_solver.Mult(z, y);
+}
+
+void IM_Evolution::ImplicitSolve(const double _dt, const Vector &x, Vector &y)
+{
+   this->initA(_dt);
+
+   // y = (M + dt S - dt K)^{-1} (-S x + K x + b)
+   S.Mult(-1.0, x, 0.0, z);
+   K.Mult(1.0, x, 1.0, z);
+   z += b;
+   A_solver->Mult(z, y);
+}
+
+// Implementation of class IMEX_Evolution
+IMEX_Evolution::IMEX_Evolution(HypreParMatrix &_M, HypreParMatrix &_S,
+                               HypreParMatrix &_K, const Vector &_b)
+   : TimeDependentOperator(_M.Height()),
+     M(_M), S(_S), K(_K), A(NULL), b(_b),
+     M_prec(M), M_solver(M.GetComm()),
+     A_prec(NULL), A_solver(NULL), dt(-1.0), z(M.Height())
+{
+   M_prec.SetType(HypreSmoother::Jacobi);
+   M_solver.SetPreconditioner(M_prec);
+   M_solver.SetOperator(M);
+
+   M_solver.iterative_mode = false;
+   M_solver.SetRelTol(1e-9);
+   M_solver.SetAbsTol(0.0);
+   M_solver.SetMaxIter(100);
+   M_solver.SetPrintLevel(0);
+}
+
+void IMEX_Evolution::initA(double _dt)
+{
+   if (fabs(dt - _dt) > 1e-4 * _dt)
+   {
+      delete A_solver;
+      delete A_prec;
+      delete A;
+
+      A = Add(_dt, S, 1.0, M); // A = M + dt * S
+      dt = _dt;
+
+      A_prec = new HypreBoomerAMG(*A);
+      A_solver = new CGSolver(A->GetComm());
+      A_solver->SetOperator(*A);
+      A_solver->SetPreconditioner(*A_prec);
+
+      A_solver->iterative_mode = false;
+      A_solver->SetRelTol(1e-9);
+      A_solver->SetAbsTol(0.0);
+      A_solver->SetMaxIter(100);
+      A_solver->SetPrintLevel(0);
+   }
+}
+
+void IMEX_Evolution::Mult(const Vector &x, Vector &y) const
+{
+   // y = M^{-1} (-S x + K x + b)
+   S.Mult(-1.0, x, 0.0, z);
+   K.Mult(1.0, x, 1.0, z);
+   z += b;
+   M_solver.Mult(z, y);
+}
+
+void IMEX_Evolution::ExplicitMult(const Vector &x, Vector &y) const
+{
+   // y = M^{-1} (K x + b)
+   K.Mult(1.0, x, 0.0, z);
+   z += b;
+   M_solver.Mult(z, y);
+}
+
+void IMEX_Evolution::ImplicitSolve(const double _dt, const Vector &x, Vector &y)
+{
+   this->initA(_dt);
+   // y = (M + dt S)^{-1} (-S x + b)
+   S.Mult(-1.0, x, 0.0, z);
+   z += b;
+   A_solver->Mult(z, y);
+}
+
+// Velocity coefficient
+void velocity_function(const Vector &x, Vector &v)
+{
+   int dim = x.Size();
+
+   // map to the reference [-1,1] domain
+   Vector X(dim);
+   for (int i = 0; i < dim; i++)
+   {
+      double center = (bb_min[i] + bb_max[i]) * 0.5;
+      X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
+   }
+
+   switch (problem)
+   {
+      case 0:
+      {
+         // Translations in 1D, 2D, and 3D
+         switch (dim)
+         {
+            case 1: v(0) = 1.0; break;
+            case 2: v(0) = sqrt(2./3.); v(1) = sqrt(1./3.); break;
+            case 3: v(0) = sqrt(3./6.); v(1) = sqrt(2./6.); v(2) = sqrt(1./6.);
+               break;
+         }
+         break;
+      }
+      case 1:
+      case 2:
+      {
+         // Clockwise rotation in 2D around the origin
+         const double w = M_PI/2;
+         switch (dim)
+         {
+            case 1: v(0) = 1.0; break;
+            case 2: v(0) = w*X(1); v(1) = -w*X(0); break;
+            case 3: v(0) = w*X(1); v(1) = -w*X(0); v(2) = 0.0; break;
+         }
+         break;
+      }
+      case 3:
+      {
+         // Clockwise twisting rotation in 2D around the origin
+         const double w = M_PI/2;
+         double d = max((X(0)+1.)*(1.-X(0)),0.) * max((X(1)+1.)*(1.-X(1)),0.);
+         d = d*d;
+         switch (dim)
+         {
+            case 1: v(0) = 1.0; break;
+            case 2: v(0) = d*w*X(1); v(1) = -d*w*X(0); break;
+            case 3: v(0) = d*w*X(1); v(1) = -d*w*X(0); v(2) = 0.0; break;
+         }
+         break;
+      }
+   }
+}
+
+// Initial condition
+double u0_function(const Vector &x)
+{
+   int dim = x.Size();
+
+   // map to the reference [-1,1] domain
+   Vector X(dim);
+   for (int i = 0; i < dim; i++)
+   {
+      double center = (bb_min[i] + bb_max[i]) * 0.5;
+      X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
+   }
+
+   switch (problem)
+   {
+      case 0:
+      case 1:
+      {
+         switch (dim)
+         {
+            case 1:
+               return exp(-40.*pow(X(0)-0.5,2));
+            case 2:
+            case 3:
+            {
+               double rx = 0.45, ry = 0.25, cx = 0., cy = -0.2, w = 10.;
+               if (dim == 3)
+               {
+                  const double s = (1. + 0.25*cos(2*M_PI*X(2)));
+                  rx *= s;
+                  ry *= s;
+               }
+               return ( erfc(w*(X(0)-cx-rx))*erfc(-w*(X(0)-cx+rx)) *
+                        erfc(w*(X(1)-cy-ry))*erfc(-w*(X(1)-cy+ry)) )/16;
+            }
+         }
+      }
+      case 2:
+      {
+         double x_ = X(0), y_ = X(1), rho, phi;
+         rho = hypot(x_, y_);
+         phi = atan2(y_, x_);
+         return pow(sin(M_PI*rho),2)*sin(3*phi);
+      }
+      case 3:
+      {
+         const double f = M_PI;
+         return sin(f*X(0))*sin(f*X(1));
+      }
+   }
+   return 0.0;
+}

examples/exSBP.cpp

@@ -1,358 +0,0 @@
-//                                MFEM Example 1
-//
-// Compile with: make exSBP
-//
-// Sample runs:  exSBP -sbp -o 0 -p 0 -r 1
-//               exSBP -sbp -o 4 -p 3
-//               
-//
-// Description:  This example code builds on Example 1 but adds SBP operators.
-//               It 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 or SBP 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;
-
-int problem;
-
-// Prescribed time-independent boundary and right-hand side functions.
-double bdr_func(const Vector &pt);
-double rhs_func(const Vector &pt);
-
-int main(int argc, char *argv[])
-{
-   // 1. Parse command-line options.
-   const char *mesh_file = "../data/unitGridTestMesh.msh";
-   int order = 1;
-   bool static_cond = false;
-	bool visualization = 1;
-	bool sbp = 0;
-	problem = 1;
-	int ref_levels = 0;
-	bool convOut = false;
-
-   OptionsParser args(argc, argv);
-   args.AddOption(&mesh_file, "-m", "--mesh",
-                  "Mesh file to use.");
-   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(&visualization, "-vis", "--visualization", "-no-vis",
-                  "--no-visualization",
-                  "Enable or disable GLVis visualization.");
-   args.AddOption(&sbp, "-sbp", "--summationbyparts", "-no-sbp",
-                  "--no-summationbyparts",
-                  "Enable or disable use of SBP operators.");
-	args.AddOption(&problem, "-p", "--problem",
-                  "Problem setup to use: 0 = transcendental manufactured solution, "
-						"1 = linear displacement, "
-						"2 = quadratic displacement, "
-						"3 = cubic displacement, "
-						"4 = quartic displacement.");
-	args.AddOption(&ref_levels, "-r", "--ref-levels",
-                  "Number of initial uniform refinement levels.");
-   args.Parse();
-   if (!args.Good())
-   {
-      args.PrintUsage(cout);
-      return 1;
-   }
-   args.PrintOptions(cout);
-
-   // 2. 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();
-
-   // 3. Refine the mesh to increase the resolution. In this example we do
-   //    'ref_levels' of uniform refinement. We choose 'ref_levels' to be the
-   //    largest number that gives a final mesh with no more than 50,000
-   //    elements.
-   {
-      for (int l = 0; l < ref_levels; l++)
-      {
-         mesh->UniformRefinement();
-      }
-   }
-
-   // 4. 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 (sbp)
-   {
-      fec = new C_SBPCollection(order, dim);
-   }
-   else if (order > 0)
-   {
-      fec = new H1_FECollection(order, dim);
-   }
-   else if (mesh->GetNodes())
-   {
-      fec = mesh->GetNodes()->OwnFEC();
-      cout << "Using isoparametric FEs: " << fec->Name() << endl;
-   }
-   else
-   {
-      fec = new H1_FECollection(order = 1, dim);
-   }
-   FiniteElementSpace *fespace = new FiniteElementSpace(mesh, fec);
-   cout << "Number of finite element unknowns: "
-        << fespace->GetTrueVSize() << endl;
-
-   // 7. 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;
-
-	// Create function coefficient bdr which holds the exact solution and is
-	// used to strongly impose boundary conditions.
-	FunctionCoefficient bdr(bdr_func);
-
-   // 5. 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;
-
-		// Project boundary conditions onto grid function to strongly impose
-		// boundary conditions. BC's are defined in the function `bdr_func`.
-		x.ProjectBdrCoefficient(bdr, ess_bdr);
-
-      fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
-   }
-
-   // 6. Set up the linear form b(.) which corresponds to the right-hand side of
-   //    the FEM linear system, which in this case is (1,phi_i) where phi_i are
-   //    the basis functions in the finite element fespace.
-   LinearForm *b = new LinearForm(fespace);
-	
-	FunctionCoefficient rhs(rhs_func);
-	b->AddDomainIntegrator(new DomainLFIntegrator(rhs));
-
-	if (problem < 0 || problem > 4)
-	{
-		mfem::out << "Invalid problem type: " << problem << "\n";
-		delete mesh;
-      return 3;
-	}
-	
-   // Start timing
-   std::chrono::time_point<std::chrono::high_resolution_clock> start = std::chrono::high_resolution_clock::now();
-
-   b->Assemble();
-
-   // End timing and compute interval
-   std::chrono::time_point<std::chrono::high_resolution_clock> finish = std::chrono::high_resolution_clock::now();
-   std::chrono::duration<double> elapsed = finish - start;
-   std::cout << "\nb->Assemble() elapsed time: " << elapsed.count() << " s\n";
-
-   // 8. 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);
-	ConstantCoefficient one(1.0);
-   a->AddDomainIntegrator(new DiffusionIntegrator(one));
-
-   // 9. 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(); }
-
-   // Start timing
-   start = std::chrono::high_resolution_clock::now();
-
-   a->Assemble();
-
-   // End timing and compute interval
-   finish = std::chrono::high_resolution_clock::now();
-   elapsed = finish - start;
-   std::cout << "\na->Assemble() elapsed time: " << elapsed.count() << " s\n";
-
-   SparseMatrix A;
-   Vector B, X;
-   a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);
-
-   mfem::out << "Size of linear system: " << A.Height() << endl;
-
-#ifndef MFEM_USE_SUITESPARSE
-   // 10. Define a simple symmetric Gauss-Seidel preconditioner and use it to
-   //     solve the system A X = B with PCG.
-   GSSmoother M(A);
-   PCG(A, M, B, X, 1, 200, 1e-12, 0.0);
-#else
-   // 10. If MFEM was compiled with SuiteSparse, use UMFPACK to solve the system.
-   UMFPackSolver umf_solver;
-   umf_solver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
-   umf_solver.SetOperator(A);
-   umf_solver.Mult(B, X);
-#endif
-
-   // 11. Recover the solution as a finite element grid function.
-   a->RecoverFEMSolution(X, *b, x);
-
-	// 12. Compute and print the L^2 norm of the error.
-   mfem::out << "\n|| u_h - u ||_{L^2} = " << x.ComputeL2Error(bdr) << '\n' << endl;
-	// mfem::out << "h: " << 0.1 / pow(2, ref_levels) << "\n";
-
-   // 12. Save the refined mesh and the solution. This output can be viewed later
-   //     using GLVis: "glvis -m refined.mesh -g sol.gf".
-   ofstream mesh_ofs("refined.mesh");
-   mesh_ofs.precision(8);
-   mesh->Print(mesh_ofs);
-   // mesh->PrintVTK(mesh_ofs); 
-   ofstream sol_ofs("sol.gf");
-   sol_ofs.precision(8);
-   x.Save(sol_ofs);
-
-   // Save solution mesh in vtk file
-   char solFileName[32];
-   if (sbp)
-   {
-      snprintf(solFileName, 32, "exSBP_SBP_O%d_P%d.vtk", order, problem);
-   }
-   else
-   {
-      snprintf(solFileName, 32, "exSBP_FE_O%d_P%d.vtk", order, problem);
-   }
-
-	if (convOut)
-	{
-		// Save convergence study information in output file
-		char outfileName[32];
-		if (problem == 0)
-		{
-			snprintf(outfileName, 32, "convOutputP%d_manufactured.txt", order);
-		}
-		else if (problem == 1)
-		{
-			snprintf(outfileName, 32, "convOutputP%d_lin.txt", order);
-		}
-		else if (problem == 2)
-		{
-			snprintf(outfileName, 32, "convOutputP%d_quad.txt", order);
-		}
-		else if (problem == 3)
-		{
-			snprintf(outfileName, 32, "convOutputP%d_cubic.txt", order);
-		}
-		else if (problem == 4)
-		{
-			snprintf(outfileName, 32, "convOutputP%d_quartic.txt", order);
-		}
-		
-
-		ofstream outputFile;
-		outputFile.open(outfileName, ios::out | ios::app); 
-
-		if (outputFile.is_open())
-		{
-			outputFile << x.ComputeL2Error(bdr) << ", " << 0.1 / pow(2, ref_levels) << "\n";
-		}
-		outputFile.close();
-	}
-
-   ofstream omesh(solFileName);
-   omesh.precision(14); 
-   mesh->PrintVTK(omesh, 1); 
-   x.SaveVTK(omesh, "sol", 1);
-
-   // 13. 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;
-   }
-
-   // 14. Free the used memory.
-   delete a;
-   delete b;
-   delete fespace;
-   if (order > 0) { delete fec; }
-   delete mesh;
-
-   return 0;
-}
-
-// Exact solution, used for the Dirichlet BC.
-double bdr_func(const Vector &pt)
-{
-	double x = pt(0), y = pt(1), z = 0.0;
-
-	if (problem == 0) // manufactured solution
-	{
-		z = sin(M_PI*x)*sin(M_PI*y);
-	}
-	else if (problem == 1) // linear displacement
-	{
-		z = 0.5*x + 0.5*y;
-	}
-	else if (problem == 2) // quadratic displacement
-	{
-		z = 0.5*x*x + 0.5*y*y;
-	}
-	else if (problem == 3) // manufactured solution
-	{
-		z = 0.5*x*x*x + 0.5*y*y*y;
-	}
-	else if (problem == 4) // manufactured solution
-	{
-		z = 0.5*x*x*x*x + 0.5*y*y*y*y;
-	}
-	return z;
-}
-
-// right hand side function for manufactured solution
-double rhs_func(const Vector &pt)
-{
-	double x = pt(0), y = pt(1), z = 0.0;
-	if (problem == 0)
-	{
-		z = 2*M_PI*M_PI*sin(M_PI*x)*sin(M_PI*y);
-	}
-	else if (problem == 1)
-	{
-		z = 0;
-	}
-	else if (problem == 2)
-	{
-		z = -2;
-	}
-	else if (problem == 3)
-	{
-		z = -3*(x+y);
-	}
-	else if (problem == 4)
-	{
-		z = -6*(x*x + y*y);
-	}
-	return z;
-}

examples/makefile

@@ -22,9 +22,9 @@ MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 
 SEQ_EXAMPLES = ex1 ex2 ex3 ex4 ex5 ex6 ex7 ex8 ex9 ex10 ex14 ex15 ex16 ex17\
-  ex18 ex19 ex20 ex22 exSBP
+  ex18 ex19 ex20 ex22 ex23
 PAR_EXAMPLES = ex1p ex2p ex3p ex4p ex5p ex6p ex7p ex8p ex9p ex10p ex11p ex12p\
- ex13p ex14p ex15p ex16p ex17p ex18p ex19p ex20p ex22p
+ ex13p ex14p ex15p ex16p ex17p ex18p ex19p ex20p ex22p ex23p
 
 ifeq ($(MFEM_USE_MPI),NO)
    EXAMPLES = $(SEQ_EXAMPLES)
@@ -117,7 +117,7 @@ clean-build:
 
 clean-exec:
 	@rm -f refined.mesh displaced.mesh mesh.* ex5.mesh
-	@rm -rf Example5* Example9* Example15* Example16*
+	@rm -rf Example5* Example9* Example15* Example16* Example23*
 	@rm -f sphere_refined.* sol.* sol_u.* sol_p.*
 	@rm -f ex9.mesh ex9-mesh.* ex9-init.* ex9-final.*
 	@rm -f deformed.* velocity.* elastic_energy.* mode_*
@@ -126,3 +126,4 @@ clean-exec:
 	@rm -f deformation.* pressure.*
 	@rm -f ex20.dat ex20p_?????.dat gnuplot_ex20.inp gnuplot_ex20p.inp
 	@rm -f ex22*.mesh ex22*.sol ex22p_*.*
+	@rm -f ex23.mesh ex23-mesh.* ex23-init.* ex23-final.*

fem/bilininteg.cpp

@@ -417,8 +417,7 @@ void DiffusionIntegrator::AssembleElementMatrix
       }
       else
       {
-         // ir = &IntRules.Get(el.GetGeomType(), order);
-         ir = &IntRules.Get(el, order); // New overload
+         ir = &IntRules.Get(el.GetGeomType(), order);
       }
    }
 

fem/fe.cpp

@@ -31,7 +31,6 @@ FiniteElement::FiniteElement(int D, Geometry::Type G, int Do, int O, int F)
    DerivType = NONE;
    DerivRangeType = SCALAR;
    DerivMapType = VALUE;
-   OperatorType = FE;
    for (int i = 0; i < Geometry::MaxDim; i++) { Orders[i] = -1; }
 #ifndef MFEM_THREAD_SAFE
    vshape.SetSize(Dof, Dim);
@@ -11846,480 +11845,6 @@ void NURBS3DFiniteElement::CalcDShape(const IntegrationPoint &ip,
    }
 }
 
-/// SBP_SegmentElement is a segment element with nodes at Gauss Lobatto
-/// points with ordering consistent with SBP_TriangleElement's edges.
-
-//////////////////////////////////////////////////////////////////////////
-/// Not currently implemented as collocated SBP type element
-//////////////////////////////////////////////////////////////////////////
-SBP_SegmentElement::SBP_SegmentElement(const int p)
-   : NodalTensorFiniteElement(1, p+1, BasisType::GaussLobatto, H1_DOF_MAP)
-{
-   const double *cp = poly1d.ClosedPoints(p+1, b_type);
-
-#ifndef MFEM_THREAD_SAFE
-   shape_x.SetSize(p+2);
-   dshape_x.SetSize(p+2);
-#endif
-
-   Nodes.IntPoint(0).x = cp[0];
-   Nodes.IntPoint(1).x = cp[p+1];
-
-   switch (p)
-   {
-      case 1:
-         Nodes.IntPoint(2).x = cp[1];
-         break;
-      case 2:
-         Nodes.IntPoint(2).x = cp[1];
-         Nodes.IntPoint(3).x = cp[2];
-         break;
-      case 3:
-         Nodes.IntPoint(2).x = cp[2];
-         Nodes.IntPoint(3).x = cp[1];
-         Nodes.IntPoint(4).x = cp[3];
-         break;
-      case 4:
-         Nodes.IntPoint(2).x = cp[2];
-         Nodes.IntPoint(3).x = cp[3];
-         Nodes.IntPoint(4).x = cp[1];
-         Nodes.IntPoint(5).x = cp[4];
-         break;
-   }
-}
-
-void SBP_SegmentElement::CalcShape(const IntegrationPoint &ip,
-                                  Vector &shape) const
-{
-   const int p = Order;
-
-#ifdef MFEM_THREAD_SAFE
-   Vector shape_x(p+2);
-#endif
-
-   basis1d.Eval(ip.x, shape_x);
-
-   shape(0) = shape_x(0);
-   shape(1) = shape_x(p+1);
-
-   switch (p)
-   {
-      case 1:
-         shape(2) = shape_x(1);
-         break;
-      case 2:
-         shape(2) = shape_x(1);
-         shape(3) = shape_x(2);
-         break;
-      case 3:
-         shape(2) = shape_x(2);
-         shape(3) = shape_x(1);
-         shape(4) = shape_x(3);
-         break;
-      case 4:
-         shape(2) = shape_x(2);
-         shape(3) = shape_x(3);
-         shape(4) = shape_x(1);
-         shape(5) = shape_x(4);
-         break;
-   }
-}
-
-void SBP_SegmentElement::CalcDShape(const IntegrationPoint &ip,
-                                   DenseMatrix &dshape) const
-{
-   const int p = Order;
-
-#ifdef MFEM_THREAD_SAFE
-   Vector shape_x(p+2), dshape_x(p+2);
-#endif
-
-   basis1d.Eval(ip.x, shape_x, dshape_x);
-
-   dshape(0,0) = dshape_x(0);
-   dshape(1,0) = dshape_x(p+1);
-
-   switch (p)
-   {
-      case 1:
-         dshape(2,0) = dshape_x(1);
-         break;
-      case 2:
-         dshape(2,0) = dshape_x(1);
-         dshape(3,0) = dshape_x(2);
-         break;
-      case 3:
-         dshape(2,0) = dshape_x(2);
-         dshape(3,0) = dshape_x(1);
-         dshape(4,0) = dshape_x(3);
-         break;
-      case 4:
-         dshape(2,0) = dshape_x(2);
-         dshape(3,0) = dshape_x(3);
-         dshape(4,0) = dshape_x(1);
-         dshape(5,0) = dshape_x(4);
-         break;
-   }
-}
-
-// Leftover function from H1_Segment element
-// void SBP_SegmentElement::ProjectDelta(int vertex, Vector &dofs) const
-// {
-//    const int p = Order;
-//    const double *cp = poly1d.ClosedPoints(p, b_type);
-
-//    switch (vertex)
-//    {
-//       case 0:
-//          dofs(0) = poly1d.CalcDelta(p, (1.0 - cp[0]));
-//          dofs(1) = poly1d.CalcDelta(p, (1.0 - cp[p]));
-//          for (int i = 1; i < p; i++)
-//          {
-//             dofs(i+1) = poly1d.CalcDelta(p, (1.0 - cp[i]));
-//          }
-//          break;
-
-//       case 1:
-//          dofs(0) = poly1d.CalcDelta(p, cp[0]);
-//          dofs(1) = poly1d.CalcDelta(p, cp[p]);
-//          for (int i = 1; i < p; i++)
-//          {
-//             dofs(i+1) = poly1d.CalcDelta(p, cp[i]);
-//          }
-//          break;
-//    }
-// }
-
-SBP_TriangleElement::SBP_TriangleElement(const int p, const int Do)
-   : NodalFiniteElement(2, Geometry::TRIANGLE, Do, p,
-                        FunctionSpace::Pk)
-{
-   // Set operator type to SBP to be used in IntegrationRules::Get()
-   OperatorType = SBP;
-
-   // Create Dx and Dy matrixes
-   Dx = new DenseMatrix(Dof);
-   Dy = new DenseMatrix(Dof);
-
-   // Data to be stored in Dx and Dy matrices depending upon order
-
-   // There is probably a better way to input these and constuct Dx and Dy,
-   // I did it this way because it was simple and fast to get the code working.
-   const double p0Dx[9] = {-0.9999999999999984,-1.000000000000001,-0.9999999999999988,
-         1.000000000000001,0.9999999999999974,0.9999999999999994,
-         1.3322676295501878e-15,-1.9984014443252818e-15,9.992007221626409e-16};
-   const double p0Dy[9] = {-0.9999999999999974,-0.9999999999999994,-1.0000000000000009,
-         1.9984014443252818e-15,-9.992007221626409e-16,-1.7486012637846216e-15,
-         1.0000000000000009,0.9999999999999991,0.9999999999999984};
-   const double p1Dx[49] = {-3.333333333333333, 0.21647921352995003, 0.10823960676497299, -0.8824293926518367, 0.7863909022744312, -0.051249687578105226, -0.362809297675581,
--0.21647921352995003, 3.333333333333333, -0.10823960676497392, 0.8824293926518368, 0.05124968757810527, -0.7863909022744306, 0.3628092976755813,
--0.10823960676497299, 0.10823960676497392, 0.0, 7.406681257404114e-16, -0.8311797050737315, 0.8311797050737307, -1.9769834191462048e-16,
-2.3531450470715645, -2.353145047071565, -1.9751150019744302e-15, 0.0, -0.2351412146963257, 0.23514121469632399, 3.1540140465066577e-17,
--2.09704240606515, -0.13666583354161405, 2.2164792135299507, 0.2351412146963257, 5.0, 0.47028242939265, 1.2743814046488378,
-0.13666583354161393, 2.0970424060651482, -2.2164792135299485, -0.23514121469632399, -0.47028242939265, -5.0, -1.2743814046488378,
-3.2652836790802304, -3.2652836790802326, 1.7792850772315852e-15, -1.0644797406959973e-16, -4.301037240689829, 4.301037240689829, 0.0};
-   const double p1Dy[49] = {-3.333333333333333, 0.10823960676497733, 0.21647921352994814, -0.0512496875781027, 0.7863909022744318, -0.8824293926518388, -0.3628092976755811,
--0.10823960676497733, 0.0, 0.10823960676497837, 0.8311797050737281, -0.8311797050737287, 9.745633233426465e-17, -3.4329014409980615e-16,
--0.21647921352994814, -0.10823960676497837, 3.333333333333333, -0.7863909022744319, 0.051249687578102346, 0.8824293926518396, 0.3628092976755813,
-0.1366658335416072, -2.216479213529942, 2.0970424060651514, -5.0, -0.4702824293926414, -0.23514121469633495, -1.274381404648838,
--2.0970424060651514, 2.216479213529943, -0.13666583354160625, 0.4702824293926414, 5.0, 0.235141214696336, 1.2743814046488384,
-2.3531450470715702, -2.5988355289137243e-16, -2.3531450470715725, 0.23514121469633495, -0.235141214696336, 0.0, 3.7848168558079887e-16,
-3.2652836790802313, 3.0896112968982564e-15, -3.2652836790802326, 4.30103724068983, -4.301037240689832, -1.2773756888351966e-15, 0.0};
-   
-   const double p2Dx[144] = {-6.654819485608504, -0.2644800033235228, -0.13224000166175998, -1.4588736121976598, 0.16568196946098376, -0.18364010800322192, -0.49629408766969996, 0.18676496181247967, -0.32850861828832895, 0.17877251381228318, -0.18297783134208934, 0.16855922332582518,
-0.2644800033235228, 6.654819485608504, 0.13224000166176486, -0.1656819694609847, 1.4588736121976598, 0.3285086182883251, -0.18676496181248323, 0.4962940876696997, 0.18364010800322153, -0.1787725138122836, -0.16855922332582482, 0.18297783134208834,
-0.13224000166175998, -0.13224000166176486, 0.0, -0.31265397966647923, 0.3126539796664765, -0.021082992351498762, -1.1303649939093317, 1.1303649939093305, 0.0210829923514981, -3.1850260779507576e-16, 0.010213290486458549, -0.010213290486458991,
-6.25007885298757, 0.7098115731146227, 1.3394662911697857, 0.0, -1.0987526276656132, 0.8391459038761706, -0.7493297423228898, -1.0608454694676046, 0.7982718911706455, 0.47711018331982885, -0.20709875867432972, -0.4372738402113571,
--0.7098115731146187, -6.25007885298757, -1.339466291169774, 1.0987526276656132, 0.0, -0.7982718911706481, 1.0608454694676077, 0.7493297423228901, -0.8391459038761684, -0.47711018331982796, 0.4372738402113582, 0.20709875867432867,
-0.7867474920341369, -1.407390435347625, 0.09032335875573284, -0.8391459038761706, 0.7982718911706481, 7.766734444360106, 0.549376313832808, -1.8999913733437686, 1.4986594846457786, 0.7173852412251129, -0.7042577654543215, 1.9425023749691277,
-2.126213783203917, 0.8001349318703544, 4.842688417639932, 0.7493297423228898, -1.0608454694676077, -0.549376313832808, 7.766734444360106, 1.5965437823412956, -1.8999913733437699, 1.7354036162948008, -0.2271475821344933, 0.2801114010137556,
--0.8001349318703392, -2.1262137832039163, -4.842688417639927, 1.0608454694676046, -0.7493297423228901, 1.8999913733437686, -1.5965437823412956, -7.766734444360106, 0.5493763138328077, -1.7354036162947999, -0.28011140101375603, 0.22714758213449346,
-1.4073904353476412, -0.7867474920341352, -0.09032335875573, -0.7982718911706455, 0.8391459038761684, -1.4986594846457786, 1.8999913733437699, -0.5493763138328077, -7.766734444360106, -0.7173852412251128, -1.9425023749691277, 0.7042577654543228,
--3.0482350464695736, 3.0482350464695807, 5.4307610872012924e-15, -1.8988845609884182, 1.8988845609884146, -2.8551722567827746, -6.906855584460779, 6.906855584460776, 2.855172256782774, 0.0, 1.0883892687387082, -1.088389268738702,
-3.11993956078517, 2.8740890922806273, -0.1741459541889496, 0.8242469961764219, -1.7403370817278376, 2.8029252875555524, 0.904040726567137, 1.114835175054942, 7.73110258063719, -1.0883892687387082, 0.0, -2.176778537477416,
--2.8740890922806335, -3.119939560785153, 0.17414595418895715, 1.7403370817278334, -0.8242469961764177, -7.73110258063719, -1.1148351750549401, -0.9040407265671375, -2.8029252875555577, 1.088389268738702, 2.176778537477416, 0.0};
-
-   const double p2Dy[144] = {-6.654819485608504, -0.13224000166176303, -0.26448000332352145, -0.3285086182883249, 0.1867649618124839, -0.49629408766970096, -0.18364010800322322, 0.16568196946098196, -1.4588736121976569, 0.16855922332582438, -0.1829778313420883, 0.17877251381228437,
-0.13224000166176303, 0.0, -0.13224000166176358, 0.021082992351497333, 1.1303649939093325, -1.13036499390933, -0.0210829923514987, 0.31265397966647596, -0.3126539796664774, -0.010213290486459194, 0.010213290486458717, 1.7550143694830702e-16,
-0.26448000332352145, 0.13224000166176358, 6.654819485608504, 0.18364010800321912, 0.49629408766969646, -0.18676496181248112, 0.32850861828832506, 1.4588736121976555, -0.16568196946098368, 0.18297783134208803, -0.16855922332582496, -0.17877251381228357,
-1.4073904353476239, -0.09032335875572672, -0.7867474920341248, -7.766734444360106, -0.5493763138328086, 1.8999913733437714, -1.498659484645771, 0.8391459038761694, -0.7982718911706419, 0.7042577654543205, -1.9425023749691286, -0.7173852412251152,
--0.8001349318703573, -4.842688417639936, -2.1262137832039025, 0.5493763138328086, -7.766734444360106, -1.5965437823412931, 1.8999913733437763, -0.7493297423228881, 1.060845469467605, 0.22714758213449282, -0.28011140101375637, -1.7354036162948017,
-2.1262137832039216, 4.842688417639926, 0.8001349318703455, -1.8999913733437714, 1.5965437823412931, 7.766734444360106, -0.5493763138328039, -1.0608454694675977, 0.7493297423228871, 0.280111401013755, -0.22714758213449193, 1.7354036162947997,
-0.7867474920341424, 0.09032335875573258, -1.4073904353476248, 1.498659484645771, -1.8999913733437763, 0.5493763138328039, 7.766734444360106, 0.7982718911706495, -0.8391459038761694, 1.942502374969129, -0.7042577654543208, 0.7173852412251135,
--0.709811573114611, -1.3394662911697717, -6.250078852987552, -0.8391459038761694, 0.7493297423228881, 1.0608454694675977, -0.7982718911706495, 0.0, 1.0987526276656179, 0.20709875867432867, 0.43727384021135907, -0.47711018331982974,
-6.250078852987558, 1.339466291169778, 0.7098115731146183, 0.7982718911706419, -1.060845469467605, -0.7493297423228871, 0.8391459038761694, -1.0987526276656179, 0.0, -0.4372738402113571, -0.20709875867432725, 0.4771101833198306,
--2.8740890922806197, 0.17414595418896062, -3.1199395607851477, -2.8029252875555484, -0.904040726567135, -1.114835175054938, -7.731102580637195, -0.8242469961764177, 1.7403370817278334, 0.0, 2.1767785374774147, 1.0883892687387102,
-3.119939560785152, -0.17414595418895248, 2.8740890922806295, 7.731102580637194, 1.1148351750549435, 0.9040407265671315, 2.8029252875555497, -1.7403370817278412, 0.8242469961764122, -2.1767785374774147, 0.0, -1.0883892687387073,
--3.048235046469594, -2.9924601909068347e-15, 3.04823504646958, 2.8551722567827835, 6.906855584460783, -6.906855584460775, -2.855172256782777, 1.8988845609884217, -1.8988845609884253, -1.0883892687387102, 1.0883892687387073, 0.0};
-
-   const double p3Dx[324] = {-10.952585131486583, -0.1691685246371182, -0.08458426231856121, 0.2782043329268357, -2.305479077267717, -0.07465619884488185, 0.08251879096720546, 0.22849304912079002, 0.5302407771816927, 0.28327860677318784, -0.11951297527437506, -0.5003055437372067, 0.12448476066856992, -0.09239160861978224, 0.10362923976480765, -0.16947318256573993, -0.12955658105892584, -0.0720914264099959,
-0.1691685246371182, 10.952585131486583, 0.08458426231855752, -0.2782043329268405, 0.0746561988448885, 2.3054790772677194, -0.2832786067731876, 0.5003055437372037, 0.11951297527437463, -0.08251879096721035, -0.5302407771816938, -0.2284930491207896, -0.12448476066856932, -0.10362923976480676, 0.09239160861978266, 0.16947318256574107, 0.07209142640999573, 0.12955658105892595,
-0.08458426231856121, -0.08458426231855752, 0.0, 1.1683416656308464e-15, 0.3017477280608998, -0.3017477280609013, -0.005074273846349578, 0.04485677642949094, -1.8051735335305155, 0.0050742738463512614, 1.805173533530519, -0.04485677642949147, 4.205982075638532e-16, 0.020855520903763268, -0.020855520903761818, 7.395777588076851e-16, -0.0973817561557449, 0.09738175615574554,
--1.8896690623176595, 1.8896690623176922, -7.935818527815355e-15, 0.0, 2.0050512471695137, -2.0050512471695145, -0.8926105190374908, 1.6696076717446249, 0.020736328416563414, 0.8926105190375015, -0.020736328416544252, -1.669607671744625, -4.429107086712402e-16, 0.4633133883374836, -0.46331338833748054, 1.783666940029547e-15, 0.7834980551346994, -0.783498055134701,
-11.456452196936608, -0.3709837064689789, -1.4994533917700335, -1.4668732071255282, 0.0, 0.43808614602067847, -0.03454227374421423, -0.36820417663301963, 0.6765153661826107, -1.2393619133905875, 1.2485385227460366, 1.2932580655616601, -0.24307113935772914, 0.052425055137118155, 0.44269118420721065, -0.48976278487794334, -0.3768721550070895, 0.17743932374896715,
-0.3709837064689459, -11.45645219693662, 1.4994533917700408, 1.4668732071255288, -0.43808614602067847, 0.0, 1.2393619133905824, -1.2932580655616563, -1.2485385227460208, 0.034542273744204106, -0.6765153661826184, 0.36820417663302096, 0.24307113935772948, -0.44269118420720915, -0.052425055137118384, 0.4897627848779442, -0.1774393237489731, 0.3768721550070885,
--0.5604988416610844, 1.924135449667926, 0.03446638735025072, 0.8926105190374908, 0.04721541624352583, -1.6940688111740314, 11.466536000042892, -1.1244968786658345, 0.8805543685036823, 1.785221038074984, -3.363676482918665, 0.026479087826965113, -0.42363907421033664, 2.460011157605557, -0.8869524625478119, 1.3603474091530658, 0.13068371284607608, 0.5768493540183615,
--1.135434157870106, -2.486132536271146, -0.22290356912810128, -1.2214664156593342, 0.36820417663301963, 1.2932580655616563, 0.8226694181207354, 12.000000000000005, -0.21904307301033732, 0.01937180635024931, 1.6167426993790468, -1.3530307323652335, -0.6132774672199158, 0.43432201350751304, -0.11093602189671813, -0.6100947628964166, 0.7813784140227625, 2.7998077074825853,
--2.634887549640154, -0.5938872755970562, 8.97031966066547, -0.01517046739396432, -0.6765153661826107, 1.2485385227460208, -0.6442037890047951, 0.21904307301033732, 12.000000000000005, -2.460828329049918, 2.586516131123318, 1.6167426993790477, -0.0585109667596005, 0.19125087414978273, -0.17058628301270465, 2.4229355524754963, 0.29161562914482175, -0.43265543914744503,
--1.9241354496679275, 0.5604988416611175, -0.034466387350262154, -0.8926105190375015, 1.6940688111740383, -0.047215416243512, -1.785221038074984, -0.026479087826961595, 3.3636764829186583, -11.466536000042892, -0.8805543685036843, 1.1244968786658245, 0.42363907421033115, 0.8869524625478146, -2.4600111576055554, -1.3603474091530607, -0.5768493540183651, -0.13068371284607755,
-0.5938872755970585, 2.634887549640159, -8.97031966066549, 0.0151704673939503, -1.2485385227460366, 0.6765153661826184, 2.460828329049923, -1.6167426993790468, -2.586516131123318, 0.6442037890047966, -12.000000000000005, -0.21904307301033812, 0.05851096675959981, 0.17058628301270434, -0.19125087414978267, -2.422935552475496, 0.43265543914744814, -0.29161562914482114,
-2.486132536271161, 1.1354341578701037, 0.22290356912810394, 1.2214664156593344, -1.2932580655616601, -0.36820417663302096, -0.019371806350251885, 1.3530307323652335, -1.6167426993790477, -0.822669418120728, 0.21904307301033812, -12.000000000000005, 0.613277467219916, 0.11093602189671813, -0.43432201350751237, 0.6100947628964153, -2.7998077074825867, -0.7813784140227636,
--4.072657651361357, 4.072657651361338, -1.3760339008438063e-14, 2.133320141015979e-15, 1.6003174861462979, -1.6003174861463, 2.0404965421712293, 4.037660156795314, 0.38522106525739636, -2.0404965421712027, -0.38522106525739186, -4.037660156795315, 0.0, -1.4579314260766258, 1.4579314260766234, 2.0161871961592337e-15, 3.273923136675281, -3.2739231366752755,
-3.0226944225627097, 3.3903460469074083, -0.6823116044539523, -2.2315915230588144, -0.34515299788282433, 2.914564209151067, -11.84886986676902, -2.859463754812241, -1.2591462686659354, -4.272088065230016, -1.1230959476442317, -0.7303740631402172, 1.4579314260766258, 0.0, 2.9158628521532504, 1.6703256795818853, -0.3883097004508776, -1.603597457093401,
--3.390346046907437, -3.022694422562723, 0.6823116044539049, 2.2315915230587997, -2.914564209151077, 0.34515299788282583, 4.272088065230004, 0.7303740631402172, 1.1230959476442337, 11.84886986676901, 1.2591462686659352, 2.8594637548122366, -1.4579314260766234, -2.9158628521532504, 0.0, -1.6703256795818826, 1.6035974570934062, 0.3883097004508795,
-3.825058001824472, -3.8250580018244977, -1.669248067139698e-14, -5.926918282616925e-15, 2.2245074239726854, -2.2245074239726894, -4.520276599333233, 2.771056460992968, -11.004997298094699, 4.520276599333215, 11.004997298094697, -2.7710564609929618, -1.3909328887646931e-15, -1.1523289737701476, 1.1523289737701459, 0.0, -0.9410755641839665, 0.9410755641839751,
-2.9241289363066816, -1.6271240279890726, 2.1979339738354153, -2.6034731278339045, 1.7117570640055604, 0.8059311674378296, -0.4342468145397358, -3.5490285023571575, -1.3245210783844836, 1.9168034714993238, -1.9651252935551593, 12.716754362100263, -2.258623293890646, 0.26788818737287284, -1.1062943201205055, 0.9410755641839665, 0.0, 1.8821511283679513,
-1.6271240279890764, -2.9241289363066834, -2.1979339738354295, 2.60347312783391, -0.8059311674378026, -1.711757064005556, -1.9168034714993114, -12.716754362100257, 1.9651252935551453, 0.4342468145397406, 1.324521078384481, 3.5490285023571624, 2.258623293890642, 1.106294320120502, -0.2678881873728741, -0.9410755641839751, -1.8821511283679513, 0.0};
-
-   const double p3Dy[324] = {-10.952585131486583, -0.08458426231856092, -0.16916852463712426, 0.28327860677318756, -0.5003055437372038, -0.11951297527437327, 0.0825187909672079, 0.5302407771816892, 0.22849304912079044, 0.2782043329268377, -0.07465619884488649, -2.305479077267721, 0.1036292397648074, -0.09239160861978217, 0.12448476066856967, -0.07209142640999439, -0.12955658105892628, -0.16947318256574065,
-0.08458426231856092, 0.0, -0.08458426231857115, 0.00507427384635043, -0.044856776429489494, 1.8051735335305168, -0.005074273846348987, -1.805173533530519, 0.04485677642949065, -2.536797215167323e-16, -0.30174772806090094, 0.3017477280609018, -0.020855520903761887, 0.020855520903762893, 3.895409306444457e-16, 0.09738175615574617, -0.09738175615574511, -2.518495690291951e-17,
-0.16916852463712426, 0.08458426231857115, 10.952585131486583, -0.08251879096720664, -0.2284930491207953, -0.5302407771816934, -0.2832786067731896, 0.11951297527436967, 0.5003055437372073, -0.2782043329268386, 2.305479077267721, 0.07465619884488589, 0.09239160861978252, -0.10362923976480695, -0.12448476066856902, 0.12955658105892542, 0.07209142640999575, 0.1694731825657405,
--1.9241354496679255, -0.034466387350256505, 0.5604988416610924, -11.466536000042892, 1.124496878665824, -0.8805543685036857, -1.7852210380749862, 3.3636764829186614, -0.026479087826939977, -0.8926105190374863, -0.04721541624350879, 1.6940688111740283, -2.4600111576055577, 0.8869524625478143, 0.42363907421033403, -0.13068371284608143, -0.5768493540183609, -1.3603474091530663,
-2.4861325362711466, 0.2229035691280941, 1.1354341578701321, -0.8226694181207276, -12.000000000000005, 0.21904307301034107, -0.01937180635025535, -1.6167426993790432, 1.3530307323652313, 1.2214664156593333, -0.3682041766330094, -1.2932580655616668, -0.4343220135075129, 0.11093602189671664, 0.6132774672199158, -0.7813784140227632, -2.799807707482588, 0.6100947628964168,
-0.5938872755970495, -8.97031966066548, 2.634887549640157, 0.6442037890047976, -0.21904307301034107, -12.000000000000005, 2.4608283290499164, -2.586516131123313, -1.6167426993790395, 0.015170467393960869, 0.6765153661826211, -1.2485385227460257, -0.19125087414978525, 0.17058628301270676, 0.058510966759598895, -0.2916156291448204, 0.43265543914744725, -2.4229355524754985,
--0.560498841661101, 0.03446638735024671, 1.9241354496679397, 1.7852210380749862, 0.026479087826969852, -3.363676482918656, 11.466536000042892, 0.8805543685036855, -1.1244968786658387, 0.8926105190374908, -1.6940688111740407, 0.047215416243521366, -0.8869524625478177, 2.460011157605558, -0.4236390742103285, 0.5768493540183687, 0.13068371284607674, 1.360347409153062,
--2.634887549640136, 8.97031966066549, -0.5938872755970316, -2.4608283290499204, 1.6167426993790432, 2.586516131123313, -0.6442037890047975, 12.000000000000005, 0.2190430730103402, -0.015170467393955334, 1.2485385227460235, -0.6765153661826085, -0.17058628301270545, 0.19125087414978356, -0.05851096675960009, -0.432655439147446, 0.2916156291448159, 2.4229355524755003,
--1.135434157870108, -0.22290356912809986, -2.4861325362711635, 0.019371806350233497, -1.3530307323652313, 1.6167426993790395, 0.8226694181207385, -0.2190430730103402, 12.000000000000005, -1.221466415659336, 1.2932580655616515, 0.3682041766330148, -0.11093602189671982, 0.43432201350751376, -0.6132774672199166, 2.7998077074825862, 0.781378414022764, -0.6100947628964163,
--1.889669062317673, 1.72308862498413e-15, 1.8896690623176793, 0.8926105190374863, -1.6696076717446235, -0.0207363284165587, -0.8926105190374908, 0.020736328416551132, 1.669607671744627, 0.0, -2.005051247169507, 2.005051247169518, -0.46331338833748403, 0.4633133883374836, -7.923728999874815e-16, -0.7834980551347012, 0.7834980551346985, 2.487535508292328e-15,
-0.37098370646896894, 1.499453391770039, -11.456452196936628, 0.03454227374420177, 0.3682041766330094, -0.6765153661826211, 1.2393619133905893, -1.2485385227460235, -1.2932580655616515, 1.4668732071255233, 0.0, -0.43808614602068285, -0.0524250551371181, -0.44269118420721065, 0.24307113935773086, 0.37687215500708976, -0.17743932374896845, 0.48976278487794267,
-11.456452196936628, -1.4994533917700434, -0.37098370646896595, -1.2393619133905802, 1.2932580655616668, 1.2485385227460257, -0.03454227374421096, 0.6765153661826085, -0.3682041766330148, -1.4668732071255315, 0.43808614602068285, 0.0, 0.44269118420721193, 0.05242505513711802, -0.24307113935772984, 0.17743932374897087, -0.37687215500708626, -0.48976278487794556,
--3.3903460469074287, 0.682311604453907, -3.022694422562719, 11.848869866769022, 2.85946375481224, 1.2591462686659523, 4.272088065230031, 1.123095947644239, 0.7303740631402283, 2.2315915230588166, 0.34515299788282394, -2.9145642091510853, 0.0, -2.9158628521532504, -1.4579314260766258, 0.3883097004508778, 1.6035974570934002, -1.6703256795818757,
-3.0226944225627075, -0.68231160445394, 3.3903460469074145, -4.272088065230015, -0.7303740631402075, -1.1230959476442475, -11.848869866769023, -1.259146268665941, -2.859463754812246, -2.2315915230588144, 2.914564209151077, -0.3451529978828235, 2.9158628521532504, 0.0, 1.4579314260766267, -1.6035974570933964, -0.3883097004508765, 1.6703256795818802,
--4.072657651361349, -1.2744265588712169e-14, 4.072657651361328, -2.040496542171217, -4.037660156795314, -0.3852210652573858, 2.04049654217119, 0.3852210652573937, 4.037660156795319, 3.816536908330339e-15, -1.6003174861463092, 1.6003174861463025, 1.4579314260766258, -1.4579314260766267, 0.0, -3.2739231366752795, 3.273923136675281, -1.5485517097979382e-15,
-1.6271240279890422, -2.1979339738354438, -2.924128936306672, 0.43424681453975356, 3.5490285023571606, 1.3245210783844776, -1.9168034714993356, 1.9651252935551495, -12.71675436210026, 2.6034731278339103, -1.7117570640055617, -0.8059311674378195, -0.26788818737287295, 1.1062943201204989, 2.2586232938906448, 0.0, -1.8821511283679437, -0.941075564183975,
-2.9241289363066914, 2.1979339738354198, -1.627124027989073, 1.9168034714993096, 12.71675436210027, -1.9651252935551553, -0.434246814539738, -1.324521078384457, -3.549028502357164, -2.6034731278339014, 0.8059311674378086, 1.7117570640055457, -1.1062943201205013, 0.26788818737287207, -2.258623293890646, 1.8821511283679437, 0.0, 0.9410755641839769,
-3.825058001824488, 5.68431650770163e-16, -3.825058001824485, 4.520276599333235, -2.7710564609929693, 11.00499729809471, -4.5202765993332195, -11.004997298094716, 2.7710564609929667, -8.265791864994906e-15, -2.2245074239726823, 2.2245074239726956, 1.152328973770141, -1.152328973770144, 1.0683192052870465e-15, 0.941075564183975, -0.9410755641839769, 0.0};
-
-   const double p4Dx[729] = {-15.28499617463146, 0.1965429454017319, 0.09827147270087869, 0.4207579292099769, -0.03011871770800639, -3.318562255325795, -0.12660138709292057, -0.20489551789966803, -0.1525100012592784, -0.02465488205886147, -0.44232233201277743, -0.16408537644014787, 0.28709832395121765, 0.06989189152237547, -0.445234967991587, 0.1446267442307163, -0.13342014677484493, 0.04884580398085625, -0.06120785634410805, -0.0065053622920859784, -0.055536979099377465, -0.05310883929020609, 0.137613494649232, 0.09694858537476729, 0.05875484376287877, 0.11266078234486253, -0.10900916485152536,
--0.1965429454017319, 15.28499617463146, -0.09827147270085763, 0.030118717708007003, -0.42075792920997546, 0.12660138709291616, 3.318562255325798, -0.28709832395121665, 0.16408537644014565, 0.4452349679915819, -0.06989189152237295, 0.15251000125928654, 0.20489551789967583, 0.44232233201277626, 0.024654882058862064, -0.14462674423071573, -0.048845803980857666, 0.13342014677484318, 0.06120785634410922, 0.055536979099377264, 0.006505362292085202, -0.13761349464923325, 0.0531088392902052, 0.10900916485152572, -0.11266078234486292, -0.05875484376287771, -0.09694858537476635,
--0.09827147270087869, 0.09827147270085763, 0.0, 0.05238551664039212, -0.05238551664039154, -0.41766744995391064, 0.41766744995391425, 0.13396665873213762, 0.13365960525875686, -0.19649327861528906, -2.87332728733421, -0.1336596052587585, -0.13396665873213948, 2.873327287334207, 0.19649327861529092, 7.071416687303329e-16, 0.0957809402498582, -0.09578094024985911, 5.095198309682238e-16, -0.005670877244732207, 0.005670877244731516, -0.0381937416118884, 0.03819374161188797, 0.024952712304370763, 0.055900325561320505, -0.05590032556132075, -0.02495271230436985,
--2.688383016538256, -0.19243998400255485, -0.3347110618046863, 0.0, -1.8584179757445634, 2.9254109394468375, 0.5544784209743959, -0.09619204170147806, 0.9621335517719756, -0.7640240022365051, -0.17045575971749122, -0.32599503457469536, 1.5819839652391803, -0.1673484738314025, -2.1466117678065375, 0.20084022575022525, 0.8121668775566188, 0.28698839120282377, 0.22242408399761535, 0.1007642564284906, 0.1396610379722831, -0.4047350379804057, -0.19103788701331148, -0.4861845324334061, 0.0708248987102739, -0.5390429074476013, 0.02706866117436579,
-0.1924399840025509, 2.688383016538247, 0.3347110618046826, 1.8584179757445634, 0.0, -0.5544784209744025, -2.9254109394468055, -1.5819839652391763, 0.32599503457469153, 2.1466117678065477, 0.1673484738314174, -0.9621335517719789, 0.09619204170146632, 0.1704557597174755, 0.7640240022365011, -0.20084022575022814, -0.28698839120282255, -0.8121668775566167, -0.22242408399761296, -0.13966103797228493, -0.10076425642849256, 0.19103788701331245, 0.40473503798040156, -0.027068661174367476, 0.5390429074475989, -0.07082489871027405, 0.48618453243340803,
-16.79909075460827, -0.6408763879659072, 2.114299162464476, -2.3177353997619106, 0.43930042352937404, 0.0, 0.23758827371798064, 0.14587872928125803, -0.026683834684955113, 0.27594927228434313, -0.8508769996321647, 1.733847812899358, -1.9338802842578322, -1.4125593711637086, 1.7702761450361215, 0.49385180688379604, -0.48895648790844665, -0.15872007418498266, -0.09920419307861998, 0.12317929238116793, 0.08585557349611746, 0.09370614907316197, -0.2607460957437406, -0.00372051314937977, -0.04020180712625255, 0.5265453421646704, -0.4929563256972127,
-0.6408763879659295, -16.79909075460828, -2.1142991624644942, -0.4393004235293688, 2.3177353997618853, -0.23758827371798064, 0.0, 1.9338802842578542, -1.733847812899381, -1.7702761450361162, 1.4125593711636986, 0.026683834684962007, -0.14587872928127077, 0.8508769996321606, -0.27594927228435256, -0.4938518068837918, 0.15872007418498413, 0.4889564879084464, 0.0992041930786239, -0.08585557349611386, -0.1231792923811642, 0.26074609574374336, -0.09370614907315725, 0.49295632569721404, -0.5265453421646692, 0.040201807126253306, 0.0037205131493783417,
-1.3091556741914314, 1.834380779552396, -0.8559641188321434, 0.09619204170147806, 1.5819839652391763, -0.1841259491984436, -2.440914756589403, 19.910396398833793, 0.9292089878722896, -0.6009288150226979, -0.17768908386483087, -0.22980299287322684, -1.9242671035439767, 2.9524607667234655, 0.2041356965652301, -1.1190995177139584, -0.24116380960165718, 0.7245564881296083, -0.49029775323051117, -0.04870816987800388, -0.7132041988172136, -0.454608663355292, 0.08758856822375095, 1.0907550537346067, 4.030854005025503, 0.7748473349352286, 0.07689148379455818,
-0.9744446123867609, -1.0484041028347049, -0.8540022369858385, -0.9621335517719756, -0.32599503457469153, 0.03367993684774262, 2.18843676448695, -0.9292089878722896, 19.910396398833793, 0.3767893371095702, 2.3244821244241383, 3.163967930478371, -0.22980299287321784, -4.587526524395932, -0.3514744230298633, 1.53672336568623, -0.04032358385143443, -0.832111126511136, -0.6124399423887177, 0.1737159141196091, -0.3506367152582221, 0.14771638250483257, 0.2886628025018225, 3.839816118012189, 0.6860200157542047, 0.11465722939811557, -0.9936515708028879,
-0.12480694029056734, -2.253850330037397, 0.9946802760235417, 0.6053185391404056, -1.7007108357677436, -0.27594927228434313, 1.7702761450361162, 0.47610199597409236, -0.29852147371705867, 17.141859608597226, -0.11879413685898123, -0.27846505302304264, 0.16173198913849687, -1.6885086434480685, 1.7017539992643287, 0.5994748509692391, -0.6941668393814941, 3.674640424508895, -0.09744008746280787, 0.20279671720873968, -0.3963039467276607, -0.7855375215315802, 0.6600584756782935, -0.4640248521828324, 0.7038226731025359, -0.1687199363657259, 0.20287631854527582,
-2.2391061027550703, 0.35380388805763474, 14.545240424570853, 0.13504815445355628, -0.13258632374176726, 0.8508769996321647, -1.4125593711636986, 0.14077894981231426, -1.8416334037878173, 0.11879413685898123, 17.141859608597226, -3.634591120025567, 2.339166352039784, 3.5405522900722497, -1.6885086434480694, 3.1856839366004506, -0.20031503249770058, 0.44075477678425123, -0.2731246543464956, 0.10359252413012018, -0.0115845139666891, 0.16267451141902461, -0.1724404495151067, 0.4430765773587932, -0.3703187031096739, 0.1671021499810791, -0.2589921793669097,
-1.048404102834719, -0.9744446123868129, 0.8540022369858489, 0.32599503457469536, 0.9621335517719789, -2.1884367644869207, -0.03367993684775132, 0.22980299287322684, -3.163967930478371, 0.3514744230298829, 4.587526524395912, -19.910396398833793, 0.9292089878722776, -2.324482124424116, -0.37678933710957857, -1.536723365686234, 0.8321111265111332, 0.04032358385143004, 0.612439942388719, 0.35063671525822543, -0.17371591411960696, -0.28866280250182014, -0.14771638250483188, 0.9936515708028864, -0.11465722939811639, -0.6860200157542059, -3.8398161180121906,
--1.8343807795524023, -1.3091556741914812, 0.8559641188321553, -1.5819839652391803, -0.09619204170146632, 2.440914756589375, 0.18412594919845968, 1.9242671035439767, 0.22980299287321784, -0.2041356965652155, -2.9524607667234513, -0.9292089878722776, -19.910396398833793, 0.17768908386480717, 0.6009288150227365, 1.1190995177139624, -0.7245564881296058, 0.24116380960166542, 0.49029775323050956, 0.7132041988172101, 0.04870816987800567, -0.08758856822375086, 0.4546086633552912, -0.07689148379455372, -0.7748473349352273, -4.030854005025498, -1.0907550537346087,
--0.35380388805764745, -2.2391061027550645, -14.545240424570837, 0.13258632374175544, -0.13504815445354384, 1.4125593711637086, -0.8508769996321606, -2.3391663520397947, 3.6345911200255827, 1.6885086434480685, -3.5405522900722497, 1.8416334037877997, -0.1407789498122955, -17.141859608597226, -0.11879413685897054, -3.1856839366004435, -0.44075477678425035, 0.20031503249770535, 0.2731246543464969, 0.01158451396669102, -0.10359252413011845, 0.17244044951510712, -0.16267451141902445, 0.25899217936691005, -0.16710214998107684, 0.37031870310967413, -0.44307657735879563,
-2.2538503300374226, -0.12480694029057034, -0.994680276023551, 1.7007108357677354, -0.6053185391404025, -1.7702761450361215, 0.27594927228435256, -0.16173198913850842, 0.2784650530230271, -1.7017539992643287, 1.6885086434480694, 0.29852147371706533, -0.476101995974123, 0.11879413685897054, -17.141859608597226, -0.5994748509692417, -3.6746404245088935, 0.6941668393814934, 0.09744008746280555, 0.39630394672766056, -0.2027967172087368, -0.6600584756782925, 0.7855375215315809, -0.20287631854527904, 0.16871993636572646, -0.7038226731025339, 0.4640248521828319,
--3.766522549523058, 3.7665225495230428, -1.8416130814168517e-14, -0.8186227860390011, 0.818622786039013, -2.5406965754889375, 2.5406965754889157, 4.561438634236725, -6.263669333621437, -3.0840905707320223, -16.389240973690264, 6.263669333621453, -4.561438634236742, 16.38924097369023, 3.0840905707320356, 0.0, 0.8064071592627126, -0.8064071592627154, -3.27587819047868e-15, -0.2376302079323201, 0.23763020793230094, -0.5181032452374912, 0.5181032452374866, -2.3969263418417017, 0.8425636809391238, -0.8425636809391335, 2.396926341841694,
-3.474668492754472, 1.2720940592425334, -2.494428490280513, -3.3103842098885403, 1.169761861647588, 2.5155118541145702, -0.8165598329744024, 0.9829813174648175, 0.1643585314258267, 3.571248068855471, 1.030551493366546, -3.3916767725891264, 2.953285123732874, 2.267530737757596, 18.904752827804817, -0.8064071592627126, 0.0, -1.6128143185254504, 0.281848015316382, -1.069973391943026, 0.5194782232486872, 1.842445782443736, -3.2355644271100483, -0.020449207033394267, -0.5385524522708879, -0.8386380852683435, 0.9998821015046064,
--1.2720940592424965, -3.474668492754426, 2.4944284902805367, -1.169761861647593, 3.310384209888532, 0.8165598329743948, -2.515511854114569, -2.953285123732884, 3.391676772589138, -18.904752827804824, -2.2675307377576006, -0.1643585314258088, -0.9829813174648512, -1.0305514933665705, -3.571248068855467, 0.8064071592627154, 1.6128143185254504, 0.0, -0.2818480153163751, -0.5194782232487007, 1.0699733919430288, 3.235564427110037, -1.8424457824437297, -0.9998821015046048, 0.8386380852683486, 0.5385524522708905, 0.020449207033400366,
-2.2431057877263694, -2.243105787726412, -1.867255202307378e-14, -1.275749958493393, 1.2757499584933791, 0.7181857904438594, -0.7181857904438878, 2.8121834969091495, 3.512750134965057, 0.7054146005697136, 1.977277771072542, -3.5127501349650645, -2.8121834969091406, -1.9772777710725515, -0.705414600569697, 4.609760570128499e-15, -0.39661177621033, 0.3966117762103203, 0.0, 1.473286325955928, -1.4732863259559217, 0.03294149361242836, -0.0329414936124289, -3.953179136944931, -1.408551582879458, 1.4085515828794593, 3.9531791369449287,
-0.23840429448465078, -2.035283159571775, 0.20782262815464633, -0.5779499847581375, 0.8010488801122674, -0.8917528051961954, 0.6215488579847753, 0.27937372870475113, -0.996376229788628, -1.468140772361843, -0.749954981917987, -2.0111346167968778, -4.090700119723205, -0.08386574258491672, -2.869030576268716, 0.3343892230125261, 1.5056485212427524, 0.7310009992228542, -1.473286325955928, 0.0, -2.9465726519118554, -0.8795833962368981, 0.19919319766269372, 0.49291247512221925, 0.5258539687346377, 4.152372334607633, 0.5289681866425602,
-2.035283159571782, -0.23840429448462236, -0.20782262815462102, -0.801048880112257, 0.5779499847581487, -0.6215488579848013, 0.8917528051961683, 4.090700119723225, 2.0111346167968587, 2.869030576268717, 0.0838657425849028, 0.9963762297886157, -0.27937372870476146, 0.7499549819179745, 1.468140772361822, -0.3343892230124991, -0.7310009992228351, -1.5056485212427564, 1.4732863259559217, 2.9465726519118554, 0.0, -0.19919319766269705, 0.8795833962369017, -0.5289681866425588, -4.152372334607625, -0.5258539687346414, -0.49291247512221426,
-1.759402157237792, 4.558892315982919, 1.2652912826383578, 2.0985063083979236, -0.9905102682507421, -0.6132401678100317, -1.706397937502056, 2.357095527641503, -0.7658930694187117, 5.140784956908139, -1.0645890975962542, 1.4966846336620354, 0.4541370173391847, -1.1285002237720632, 4.319613754198333, 0.6590565197748737, -2.3436948454060653, -4.115820254824023, -0.029778238966182636, 0.795120126366912, 0.1800653822873389, 0.0, 3.1163752731045724, 0.15960718397414111, -0.7346658981606436, -0.1973888096528869, -0.7116415547893532,
--4.558892315982878, -1.7594021572377623, -1.2652912826383433, 0.9905102682507371, -2.0985063083979023, 1.7063979375020377, 0.6132401678100009, -0.4541370173391852, -1.4966846336620476, -4.319613754198339, 1.1285002237720605, 0.7658930694187083, -2.357095527641499, 1.064589097596253, -5.140784956908144, -0.6590565197748679, 4.115820254824037, 2.3436948454060573, 0.029778238966183125, -0.1800653822873359, -0.7951201263669153, -3.1163752731045724, 0.0, 0.71164155478935, 0.1973888096528901, 0.7346658981606372, -0.15960718397413862,
--3.2117356080303803, -3.6112813302216713, -0.8266393399665194, 2.520812908731641, 0.14034800771819048, 0.024348115151799948, -3.2260489080343286, -5.655444047050364, -19.909020941069866, 3.0367129695880832, -2.8996213945368607, -5.151973251646951, 0.3986738203099664, -1.6949198009668083, 1.3276813620016823, 3.049025358424062, 0.026012543529633645, 1.271906370130427, 3.57356938337672, -0.4455798406197396, 0.47817325020768053, -0.15960718397414111, -0.71164155478935, 0.0, -1.5581876365522882, -0.3569959936270113, 1.46933179632127,
--1.9464443253920187, 3.7322529760164143, -1.8518791729839088, -0.3672192491087197, -2.794877724005471, 0.26309226440541583, 3.4458651559413247, -20.899531209320617, -3.556937738652459, -4.606019332038912, 2.4234728017780744, 0.594485025057365, 4.017497542393682, 1.0935648461639937, -1.1041521086202575, -1.071788475275629, 0.6850690633045105, -1.0667948963999776, 1.2732933765745933, -0.4753580795859133, 3.753634765664056, 0.7346658981606436, -0.1973888096528901, 1.5581876365522882, 0.0, -1.4232831095787026, -0.35699599362703105,
--3.7322529760164014, 1.9464443253919836, 1.8518791729839168, 2.794877724005484, 0.36721924910872045, -3.4458651559413327, -0.26309226440542083, -4.017497542393689, -0.5944850250573608, 1.1041521086202537, -1.0935648461640088, 3.5569377386524654, 20.89953120932059, -2.423472801778076, 4.6060193320388985, 1.0717884752756415, 1.0667948963999712, -0.6850690633045139, -1.2732933765745946, -3.753634765664064, 0.47535807958591664, 0.1973888096528869, -0.7346658981606372, 0.3569959936270113, 1.4232831095787026, 0.0, -1.5581876365522849,
-3.6112813302216593, 3.2117356080303487, 0.8266393399664891, -0.14034800771818173, -2.5208129087316506, 3.2260489080343198, -0.024348115151790604, -0.3986738203099895, 5.151973251646959, -1.3276813620016612, 1.694919800966806, 19.909020941069873, 5.6554440470503735, 2.8996213945368763, -3.03671296958808, -3.0490253584240525, -1.2719063701304292, -0.026012543529641403, -3.5735693833767184, -0.47817325020768164, 0.4455798406197351, 0.7116415547893532, 0.15960718397413862, -1.46933179632127, 0.35699599362703105, 1.5581876365522849, 0.0};
-
-   const double p4Dy[729] = {-15.28499617463146, 0.09827147270086124, 0.19654294540174796, 0.2870983239512125, -0.164085376440144, -0.4452349679915834, 0.06989189152237135, -0.15251000125927563, -0.2048955178996692, -0.44232233201278237, -0.024654882058863847, -0.03011871770800668, 0.42075792920998006, -0.12660138709291216, -3.3185622553257907, 0.048845803980856306, -0.13342014677484385, 0.14462674423071498, -0.05553697909937705, -0.006505362292085907, -0.06120785634410891, -0.10900916485152638, 0.11266078234486171, 0.05875484376287875, 0.09694858537476714, 0.1376134946492328, -0.05310883929020566,
--0.09827147270086124, 0.0, 0.09827147270087347, -0.13396665873213587, -0.13365960525875817, 0.19649327861528584, 2.8733272873342135, 0.1336596052587607, 0.13396665873214064, -2.873327287334207, -0.19649327861528823, -0.052385516640395494, 0.05238551664038741, 0.4176674499539139, -0.41766744995391747, -0.095780940249859, 0.09578094024985785, -9.030965889809071e-16, 0.00567087724473144, -0.0056708772447323225, -2.468104976161084e-16, -0.024952712304370517, -0.05590032556132032, 0.05590032556132042, 0.02495271230436975, 0.03819374161188785, -0.03819374161188903,
--0.19654294540174796, -0.09827147270087347, 15.28499617463146, 0.20489551789967816, 0.15251000125928516, 0.024654882058868614, 0.44232233201278004, 0.16408537644014445, -0.2870983239512186, -0.06989189152236565, 0.445234967991583, -0.42075792920997346, 0.03011871770800419, 3.3185622553257907, 0.12660138709291832, 0.13342014677484448, -0.04884580398085644, -0.14462674423071606, 0.0065053622920860695, 0.055536979099378284, 0.06120785634410899, -0.0969485853747657, -0.058754843762877665, -0.11266078234486313, 0.1090091648515255, 0.05310883929020579, -0.1376134946492322,
--1.8343807795523694, 0.8559641188321322, -1.309155674191496, -19.910396398833793, -0.9292089878722894, 0.6009288150227082, 0.1776890838648202, 0.22980299287321443, 1.9242671035439682, -2.9524607667234357, -0.20413569656523223, -0.09619204170148214, -1.581983965239188, 0.18412594919846403, 2.4409147565893874, 0.24116380960165718, -0.7245564881296103, 1.1190995177139629, 0.048708169878006975, 0.7132041988172096, 0.49029775323050717, -1.0907550537346093, -4.0308540050255, -0.774847334935228, -0.07689148379455714, 0.45460866335528777, -0.0875885682237497,
-1.0484041028346944, 0.8540022369858469, -0.9744446123868041, 0.9292089878722894, -19.910396398833793, -0.376789337109576, -2.3244821244241085, -3.163967930478372, 0.22980299287322126, 4.587526524395933, 0.35147442302985465, 0.9621335517719841, 0.32599503457468426, -0.03367993684773616, -2.1884367644869362, 0.0403235838514332, 0.8321111265111371, -1.5367233656862271, -0.17371591411960613, 0.3506367152582226, 0.6124399423887191, -3.8398161180121875, -0.686020015754205, -0.11465722939811551, 0.99365157080289, -0.14771638250483227, -0.2886628025018202,
-2.2538503300374044, -0.9946802760235254, -0.1248069402906035, -0.4761019959741005, 0.2985214737170633, -17.141859608597226, 0.1187941368589855, 0.27846505302302493, -0.1617319891384971, 1.6885086434480594, -1.7017539992643365, -0.605318539140416, 1.700710835767743, 0.2759492722843605, -1.7702761450361306, 0.6941668393814986, -3.674640424508889, -0.5994748509692356, -0.20279671720873876, 0.3963039467276609, 0.097440087462807, 0.4640248521828354, -0.703822673102532, 0.16871993636572868, -0.20287631854527582, 0.7855375215315806, -0.6600584756782898,
--0.35380388805762664, -14.545240424570872, -2.2391061027550836, -0.14077894981230582, 1.8416334037877937, -0.1187941368589855, -17.141859608597226, 3.6345911200255987, -2.339166352039789, -3.540552290072251, 1.6885086434480605, -0.135048154453527, 0.13258632374176318, -0.8508769996321724, 1.4125593711637054, 0.20031503249770397, -0.44075477678425157, -3.1856839366004435, -0.103592524130115, 0.01158451396669216, 0.2731246543464982, -0.44307657735879197, 0.37031870310967574, -0.16710214998107722, 0.2589921793669143, -0.16267451141902012, 0.17244044951510798,
-0.9744446123867433, -0.854002236985863, -1.048404102834697, -0.22980299287321443, 3.163967930478372, -0.35147442302986054, -4.587526524395952, 19.910396398833793, -0.9292089878722678, 2.324482124424122, 0.37678933710956625, -0.32599503457469653, -0.9621335517719883, 2.188436764486939, 0.03367993684775833, -0.8321111265111345, -0.0403235838514329, 1.5367233656862256, -0.3506367152582263, 0.17371591411960866, -0.612439942388722, -0.9936515708028895, 0.11465722939811439, 0.686020015754204, 3.839816118012191, 0.28866280250181964, 0.1477163825048354,
-1.3091556741914387, -0.8559641188321627, 1.8343807795524083, -1.9242671035439682, -0.22980299287322126, 0.20413569656521577, 2.952460766723458, 0.9292089878722678, 19.910396398833793, -0.17768908386481014, -0.6009288150226963, 1.5819839652391892, 0.09619204170147512, -2.440914756589397, -0.18412594919844522, 0.7245564881296075, -0.24116380960166175, -1.119099517713959, -0.7132041988172084, -0.04870816987800363, -0.4902977532305103, 0.0768914837945529, 0.774847334935229, 4.030854005025501, 1.090755053734602, 0.08758856822374894, -0.4546086633552895,
-2.2391061027550956, 14.545240424570837, 0.3538038880575977, 2.3391663520397716, -3.634591120025584, -1.6885086434480594, 3.540552290072251, -1.8416334037878042, 0.14077894981229783, 17.141859608597226, 0.1187941368589982, -0.13258632374177015, 0.13504815445355095, -1.4125593711636983, 0.8508769996321621, 0.4407547767842533, -0.200315032497699, 3.1856839366004475, -0.011584513966690685, 0.10359252413011888, -0.27312465434649635, -0.2589921793669082, 0.16710214998107692, -0.37031870310967424, 0.4430765773587957, -0.17244044951510495, 0.16267451141902212,
-0.12480694029057937, 0.9946802760235375, -2.253850330037402, 0.16173198913851014, -0.27846505302302027, 1.7017539992643365, -1.6885086434480605, -0.29852147371705556, 0.47610199597409114, -0.1187941368589982, 17.141859608597226, -1.7007108357677325, 0.6053185391404096, 1.7702761450361277, -0.27594927228436483, 3.674640424508897, -0.694166839381494, 0.5994748509692353, -0.39630394672765873, 0.2027967172087418, -0.09744008746280586, 0.2028763185452754, -0.1687199363657282, 0.7038226731025327, -0.4640248521828344, 0.6600584756782918, -0.7855375215315795,
-0.1924399840025528, 0.33471106180470783, 2.6883830165382343, 0.09619204170148214, -0.9621335517719841, 0.7640240022365183, 0.17045575971745425, 0.32599503457469653, -1.5819839652391892, 0.16734847383142104, 2.1466117678065335, 0.0, 1.8584179757445463, -2.9254109394468206, -0.5544784209744086, -0.8121668775566203, -0.286988391202826, -0.20084022575023078, -0.10076425642849045, -0.13966103797228321, -0.2224240839976154, 0.4861845324334064, -0.07082489871027524, 0.5390429074475949, -0.027068661174365027, 0.4047350379803992, 0.1910378870133102,
--2.6883830165382765, -0.33471106180465615, -0.19243998400253687, 1.581983965239188, -0.32599503457468426, -2.146611767806547, -0.16734847383141227, 0.9621335517719883, -0.09619204170147512, -0.1704557597174845, -0.7640240022365101, -1.8584179757445463, 0.0, 0.554478420974375, 2.925410939446845, 0.2869883912028199, 0.8121668775566261, 0.20084022575022775, 0.1396610379722875, 0.10076425642848952, 0.22242408399761332, 0.027068661174365592, -0.5390429074475985, 0.07082489871027872, -0.4861845324334045, -0.191037887013308, -0.40473503798040444,
-0.640876387965887, -2.114299162464493, -16.799090754608244, -0.14587872928127424, 0.026683834684949993, -0.2759492722843605, 0.8508769996321724, -1.7338478128993726, 1.9338802842578497, 1.4125593711636983, -1.7702761450361277, 2.3177353997618972, -0.43930042352935234, 0.0, -0.23758827371795765, 0.4889564879084486, 0.15872007418498785, -0.4938518068837892, -0.12317929238116428, -0.08585557349611499, 0.09920419307862376, 0.0037205131493781296, 0.04020180712624783, -0.5265453421646701, 0.49295632569721504, -0.09370614907316001, 0.26074609574373936,
-16.799090754608244, 2.1142991624645107, -0.6408763879659182, -1.933880284257842, 1.7338478128993704, 1.7702761450361306, -1.4125593711637054, -0.02668383468496756, 0.1458787292812593, -0.8508769996321621, 0.27594927228436483, 0.4393004235293789, -2.317735399761917, 0.23758827371795765, 0.0, -0.15872007418498993, -0.4889564879084488, 0.4938518068837932, 0.08585557349611653, 0.12317929238116535, -0.0992041930786226, -0.4929563256972161, 0.52654534216467, -0.04020180712625487, -0.0037205131493803054, -0.2607460957437367, 0.09370614907315984,
--1.2720940592424979, 2.494428490280534, -3.4746684927544607, -0.9829813174648175, -0.16435853142582169, -3.5712480688554935, -1.0305514933665634, 3.3916767725891317, -2.953285123732881, -2.2675307377576113, -18.904752827804835, 3.310384209888547, -1.1697618616475773, -2.5155118541145804, 0.8165598329744322, 0.0, 1.6128143185254429, 0.806407159262733, 1.0699733919430285, -0.5194782232486954, -0.2818480153163862, 0.020449207033393636, 0.538552452270888, 0.8386380852683486, -0.9998821015046098, -1.842445782443737, 3.2355644271100434,
-3.4746684927544433, -2.4944284902805043, 1.2720940592425014, 2.9532851237328925, -3.3916767725891424, 18.904752827804796, 2.2675307377576024, 0.1643585314258205, 0.9829813174648363, 1.0305514933665378, 3.5712480688554704, 1.1697618616476022, -3.3103842098885705, -0.8165598329744215, 2.5155118541145813, -1.6128143185254429, 0.0, -0.8064071592627268, 0.5194782232486914, -1.069973391943029, 0.28184801531637504, 0.9998821015046067, -0.838638085268352, -0.5385524522708939, -0.020449207033397976, -3.2355644271100372, 1.8424457824437361,
--3.7665225495230237, 2.3519395979540512e-14, 3.7665225495230517, -4.561438634236743, 6.263669333621425, 3.084090570732004, 16.38924097369023, -6.263669333621419, 4.561438634236729, -16.38924097369025, -3.0840905707320028, 0.8186227860390236, -0.8186227860390113, 2.5406965754889024, -2.540696575488923, -0.806407159262733, 0.8064071592627268, 0.0, 0.23763020793230993, -0.23763020793231787, 4.300846213266488e-15, 2.39692634184169, -0.8425636809391248, 0.8425636809391247, -2.3969263418416866, 0.5181032452374869, -0.5181032452374932,
-2.035283159571767, -0.20782262815461824, -0.2384042944846541, -0.27937372870476895, 0.9963762297886111, 1.4681407723618363, 0.7499549819179495, 2.011134616796883, 4.090700119723195, 0.08386574258491429, 2.8690305762687025, 0.5779499847581366, -0.801048880112282, 0.891752805196169, -0.6215488579847945, -1.5056485212427562, -0.7310009992228411, -0.3343892230125118, 0.0, 2.9465726519118554, 1.4732863259559286, -0.4929124751222108, -0.5258539687346461, -4.1523723346076205, -0.5289681866425575, 0.8795833962369044, -0.1991931976626978,
-0.2384042944846482, 0.2078226281546506, -2.0352831595718124, -4.090700119723202, -2.0111346167968613, -2.8690305762687185, -0.08386574258492496, -0.9963762297886256, 0.27937372870474975, -0.7499549819179776, -1.4681407723618582, 0.8010488801122575, -0.5779499847581313, 0.6215488579847835, -0.8917528051961766, 0.7310009992228468, 1.5056485212427566, 0.334389223012523, -2.9465726519118554, 0.0, -1.47328632595593, 0.5289681866425614, 4.152372334607627, 0.5258539687346433, 0.4929124751222133, 0.19919319766269775, -0.8795833962368973,
-2.243105787726401, 9.044950905679127e-15, -2.2431057877264036, -2.812183496909127, -3.512750134965065, -0.7054146005697074, -1.9772777710725609, 3.512750134965082, 2.8121834969091446, 1.9772777710725475, 0.7054146005696991, 1.2757499584933931, -1.2757499584933811, -0.7181857904438868, 0.7181857904438784, 0.39661177621033594, -0.39661177621032023, -6.052078294524533e-15, -1.4732863259559286, 1.47328632595593, 0.0, 3.9531791369449296, 1.4085515828794586, -1.4085515828794646, -3.953179136944927, -0.03294149361242823, 0.032941493612420415,
-3.6112813302216935, 0.8266393399665112, 3.2117356080303274, 5.655444047050377, 19.909020941069855, -3.0367129695881028, 2.8996213945368523, 5.151973251646967, -0.3986738203099622, 1.6949198009667963, -1.3276813620016585, -2.520812908731642, -0.1403480077181807, -0.024348115151789213, 3.226048908034342, -0.026012543529632844, -1.2719063701304294, -3.049025358424047, 0.445579840619732, -0.4781732502076828, -3.573569383376719, 0.0, 1.5581876365522862, 0.35699599362702117, -1.4693317963212784, 0.159607183974139, 0.7116415547893488,
--3.7322529760163743, 1.8518791729839026, 1.946444325391982, 20.899531209320596, 3.556937738652461, 4.606019332038886, -2.4234728017780864, -0.5944850250573547, -4.01749754239369, -1.0935648461639944, 1.1041521086202688, 0.3672192491087266, 2.7948777240054694, -0.26309226440538497, -3.4458651559413305, -0.6850690633045106, 1.066794896399982, 1.0717884752756304, 0.4753580795859208, -3.7536347656640583, -1.273293376574594, -1.5581876365522862, 0.0, 1.4232831095787029, 0.3569959936270241, -0.7346658981606364, 0.19738880965288536,
--1.9464443253920178, -1.8518791729839061, 3.7322529760164214, 4.017497542393685, 0.5944850250573606, -1.104152108620272, 1.0935648461639964, -3.5569377386524557, -20.8995312093206, 2.4234728017780762, -4.6060193320388905, -2.7948777240054508, -0.36721924910874476, 3.445865155941331, 0.26309226440543104, -1.0667948963999776, 0.685069063304518, -1.0717884752756304, 3.753634765664053, -0.47535807958591836, 1.2732933765745995, -0.35699599362702117, -1.4232831095787029, 0.0, 1.5581876365522802, -0.19738880965288783, 0.73466589816063,
--3.211735608030375, -0.8266393399664858, -3.611281330221664, 0.39867382030998416, -5.15197325164697, 1.3276813620016612, -1.6949198009668363, -19.909020941069876, -5.65544404705034, -2.8996213945368767, 3.0367129695880966, 0.14034800771817776, 2.5208129087316324, -3.226048908034335, 0.024348115151803455, 1.2719063701304334, 0.026012543529638363, 3.0490253584240428, 0.4781732502076793, -0.4455798406197342, 3.573569383376717, 1.4693317963212784, -0.3569959936270241, -1.5581876365522802, 0.0, -0.7116415547893499, -0.1596071839741321,
--4.558892315982905, -1.2652912826383393, -1.7594021572377818, -2.3570955276414813, 0.7658930694187103, -5.140784956908142, 1.0645890975962247, -1.496684633662033, -0.4541370173391748, 1.1285002237720492, -4.319613754198328, -2.09850630839789, 0.9905102682507191, 0.6132401678100189, 1.7063979375020122, 2.3436948454060667, 4.115820254824023, -0.6590565197748682, -0.7951201263669176, -0.18006538228733954, 0.029778238966182518, -0.159607183974139, 0.7346658981606364, 0.19738880965288783, 0.7116415547893499, 0.0, -3.11637527310457,
-1.7594021572377778, 1.2652912826383786, 4.5588923159828845, 0.4541370173391787, 1.4966846336620356, 4.319613754198315, -1.1285002237720687, -0.7658930694187266, 2.35709552764149, -1.064589097596238, 5.140784956908135, -0.9905102682507304, 2.098506308397917, -1.7063979375020297, -0.6132401678100178, -4.115820254824031, -2.3436948454060658, 0.6590565197748763, 0.1800653822873396, 0.7951201263669113, -0.029778238966175454, -0.7116415547893488, -0.19738880965288536, -0.73466589816063, 0.1596071839741321, 3.11637527310457, 0.0};
-
-   // Populate the Dx and Dy matrices and create the element's Nodes
-   switch (p)
-   {
-      case 0:
-         *Dx=p0Dx;
-         *Dy=p0Dy;
-
-         Nodes.IntPoint(0).Set2w(0.0, 0.0, 0.16666666666666666);
-         Nodes.IntPoint(1).Set2w(1.0, 0.0, 0.16666666666666666);
-         Nodes.IntPoint(2).Set2w(0.0, 1.0, 0.16666666666666666);
-         break;
-      case 1:
-         *Dx=p1Dx;
-         *Dy=p1Dy;
-
-         Nodes.IntPoint(0).Set2w(0.0, 0.0, 0.024999999999999998);
-         Nodes.IntPoint(1).Set2w(1.0, 0.0, 0.024999999999999998);
-         Nodes.IntPoint(2).Set2w(0.0, 1.0, 0.024999999999999998);
-         Nodes.IntPoint(3).Set2w(0.5, 0.0, 0.06666666666666667);
-         Nodes.IntPoint(4).Set2w(0.5, 0.5, 0.06666666666666667);
-         Nodes.IntPoint(5).Set2w(0.0, 0.5, 0.06666666666666667);
-         Nodes.IntPoint(6).Set2w(0.3333333333333333, 0.3333333333333333, 0.22500000000000006);
-         break;
-      case 2:
-         *Dx=p2Dx;
-         *Dy=p2Dy;
-  
-         // vertices
-         Nodes.IntPoint(0).Set2w(0.0, 0.0, 0.006261126504899741);
-         Nodes.IntPoint(1).Set2w(1.0, 0.0, 0.006261126504899741);
-         Nodes.IntPoint(2).Set2w(0.0, 1.0, 0.006261126504899741);
-
-         // edges
-         Nodes.IntPoint(3).Set2w(0.27639320225002106, 0.0, 0.026823800250389242);
-         Nodes.IntPoint(4).Set2w(0.7236067977499789, 0.0, 0.026823800250389242);
-         Nodes.IntPoint(5).Set2w(0.7236067977499789, 0.27639320225002106, 0.026823800250389242);
-         Nodes.IntPoint(6).Set2w(0.27639320225002106, 0.7236067977499789, 0.026823800250389242);
-         Nodes.IntPoint(7).Set2w(0.0, 0.7236067977499789, 0.026823800250389242);
-         Nodes.IntPoint(8).Set2w(0.0, 0.27639320225002106, 0.026823800250389242);
-
-         // interior
-         Nodes.IntPoint(9).Set2w(0.21285435711180825, 0.5742912857763836, 0.10675793966098839);
-         Nodes.IntPoint(10).Set2w(0.21285435711180825, 0.21285435711180825, 0.10675793966098839);
-         Nodes.IntPoint(11).Set2w(0.5742912857763836, 0.21285435711180825, 0.10675793966098839);
-         break;
-      case 3:
-         *Dx=p3Dx;
-         *Dy=p3Dy;
-   
-         // vertices
-         Nodes.IntPoint(0).Set2w(0.0, 0.0, 0.0022825661430496253);
-         Nodes.IntPoint(1).Set2w(1.0, 0.0, 0.0022825661430496253);
-         Nodes.IntPoint(2).Set2w(0.0, 1.0, 0.0022825661430496253);
-
-         // edges
-         Nodes.IntPoint(3).Set2w(0.5, 0.0, 0.015504052643022513);
-         Nodes.IntPoint(4).Set2w(0.17267316464601146, 0.0, 0.011342592592592586);
-         Nodes.IntPoint(5).Set2w(0.8273268353539885, 0.0, 0.011342592592592586);
-
-         Nodes.IntPoint(6).Set2w(0.5, 0.5, 0.015504052643022513);
-         Nodes.IntPoint(7).Set2w(0.8273268353539885, 0.17267316464601146, 0.011342592592592586);
-         Nodes.IntPoint(8).Set2w(0.17267316464601146, 0.8273268353539885, 0.011342592592592586);
-
-         Nodes.IntPoint(9).Set2w(0.0, 0.5, 0.015504052643022513);
-         Nodes.IntPoint(10).Set2w(0.0, 0.8273268353539885, 0.011342592592592586);
-         Nodes.IntPoint(11).Set2w(0.0, 0.17267316464601146, 0.011342592592592586);
-
-         // interior
-         Nodes.IntPoint(12).Set2w(0.4243860251718814, 0.1512279496562372, 0.07467669469983994);
-         Nodes.IntPoint(13).Set2w(0.4243860251718814, 0.4243860251718814, 0.07467669469983994);
-         Nodes.IntPoint(14).Set2w(0.1512279496562372, 0.4243860251718814, 0.07467669469983994);
-
-         Nodes.IntPoint(15).Set2w(0.14200508409677795, 0.7159898318064442, 0.051518167995569394);
-         Nodes.IntPoint(16).Set2w(0.14200508409677795, 0.14200508409677795, 0.051518167995569394);
-         Nodes.IntPoint(17).Set2w(0.7159898318064442, 0.14200508409677795, 0.051518167995569394);
-
-         break;
-      case 4:
-         *Dx=p4Dx;
-         *Dy=p4Dy; 
-
-         // vertices
-         Nodes.IntPoint(0).Set2w(0.000000000000000000,0.000000000000000000,0.001090393904993471);
-         Nodes.IntPoint(1).Set2w(1.000000000000000000,0.000000000000000000,0.001090393904993471);
-         Nodes.IntPoint(2).Set2w(0.000000000000000000,1.000000000000000000,0.001090393904993471);
-
-         // edges
-         Nodes.IntPoint(3).Set2w(0.357384241759677534,0.000000000000000000,0.006966942871463700);
-         Nodes.IntPoint(4).Set2w(0.642615758240322466,0.000000000000000000,0.006966942871463700);
-         Nodes.IntPoint(5).Set2w(0.117472338035267576,0.000000000000000000,0.005519747637357106);
-         Nodes.IntPoint(6).Set2w(0.882527661964732424,0.000000000000000000,0.005519747637357106);
-
-         Nodes.IntPoint(7).Set2w(0.642615758240322466,0.357384241759677534,0.006966942871463700);
-         Nodes.IntPoint(8).Set2w(0.357384241759677534,0.642615758240322466,0.006966942871463700);
-         Nodes.IntPoint(9).Set2w(0.882527661964732424,0.117472338035267576,0.005519747637357106);
-         Nodes.IntPoint(10).Set2w(0.117472338035267576,0.882527661964732424,0.005519747637357106);
-
-         Nodes.IntPoint(11).Set2w(0.000000000000000000,0.642615758240322466,0.006966942871463700);
-         Nodes.IntPoint(12).Set2w(0.000000000000000000,0.357384241759677534,0.006966942871463700);
-         Nodes.IntPoint(13).Set2w(0.000000000000000000,0.882527661964732424,0.005519747637357106);
-         Nodes.IntPoint(14).Set2w(0.000000000000000000,0.117472338035267576,0.005519747637357106);
-
-         // interior
-         Nodes.IntPoint(15).Set2w(0.103677508142805172,0.792644983714389628,0.028397190663911491);
-         Nodes.IntPoint(16).Set2w(0.103677508142805172,0.103677508142805172,0.028397190663911491);
-         Nodes.IntPoint(17).Set2w(0.792644983714389628,0.103677508142805172,0.028397190663911491);
-         Nodes.IntPoint(18).Set2w(0.265331380484209678,0.469337239031580644,0.039960048027851809);
-         Nodes.IntPoint(19).Set2w(0.265331380484209678,0.265331380484209678,0.039960048027851809);
-         Nodes.IntPoint(20).Set2w(0.469337239031580644,0.265331380484209678,0.039960048027851809);
-         Nodes.IntPoint(21).Set2w(0.587085567133367348,0.088273960601581103,0.036122826526134168);
-         Nodes.IntPoint(22).Set2w(0.324640472265051494,0.088273960601581103,0.036122826526134168);
-         Nodes.IntPoint(23).Set2w(0.324640472265051494,0.587085567133367348,0.036122826526134168);
-         Nodes.IntPoint(24).Set2w(0.587085567133367348,0.324640472265051494,0.036122826526134168);
-         Nodes.IntPoint(25).Set2w(0.088273960601581103,0.324640472265051494,0.036122826526134168);
-         Nodes.IntPoint(26).Set2w(0.088273960601581103,0.587085567133367348,0.036122826526134168);
-
-         break;
-      default:
-         mfem_error("SBP elements are currently only supported for 0 <= order <= 4");
-         break;
-   }
-}
-
-/// CalcShape outputs ndofx1 vector shape based on Kronecker \delta_{i, ip}
-/// where ip is the integration point CalcShape is evaluated at. 
-void SBP_TriangleElement::CalcShape(const IntegrationPoint &ip,
-                                   Vector &shape) const
-{
-   shape = 0.0;
-   shape[ip.GetIdx()] = 1.0;
-   // for (int i = 0; i < Dof; i++)
-   // {
-   //    if (ip.x == Nodes.IntPoint(i).x && ip.y == Nodes.IntPoint(i).y)
-   //    {
-   //       shape(i) = 1;
-   //    }
-   // }
-}
-
-/// CalcDShape outputs ndof x ndim DenseMatrix dshape, where the first column
-/// is the ith row of Dx, and the second column is the ith row of Dy, where i
-/// is the integration point CalcDShape is evaluated at. Since DenseMatrices 
-/// are stored a column major we should store the transpose so accessing a row
-/// is faster, but this is not done here.
-void SBP_TriangleElement::CalcDShape(const IntegrationPoint &ip,
-                                    DenseMatrix &dshape) const
-{
-   int ipIdx = ip.GetIdx();
-   dshape = 0.0;
-
-   // for (int i = 0; i < Dof; i++)
-   // {
-   //    if (ip.x == Nodes.IntPoint(i).x && ip.y == Nodes.IntPoint(i).y)
-   //    {
-   //       ipNum = i;
-   //    }
-   // }
-
-   Vector tempVec(Dof);
-
-   // when we switch to storing Dx and Dy transpose so that access to the row we want
-   // is faster Dx->GetRow() will be replaced with Dx->GetColumnReference() or 
-   // Dx->GetColumn(), whichever is faster
-   Dx->GetRow(ipIdx, tempVec);
-   dshape.SetCol(0, tempVec);
-   Dy->GetRow(ipIdx, tempVec);
-   dshape.SetCol(1, tempVec);
-}
-
-SBP_TriangleElement::~SBP_TriangleElement()
-{
-   delete Dx;
-   delete Dy;
-}
 
 // Global object definitions
 

fem/fe.hpp

@@ -154,7 +154,6 @@ protected:
 #endif
 
 public:
-   int OperatorType;
    /// Enumeration for RangeType and DerivRangeType
    enum { SCALAR, VECTOR };
 
@@ -199,14 +198,6 @@ public:
           CURL  ///< Implements CalcCurlShape methods
         };
 
-   /** @brief Enumeration for Operator Type: defines which integration rules
-       should be called.
-    
-   */
-   enum { FE = 0, ///< Traditional Finite Element
-          SBP = 1 ///< Summation-by-parts Element with collocated integration rule
-        };
-
    /** Construct FiniteElement with given
        @param D    Reference space dimension
        @param G    Geometry type (of type Geometry::Type)
@@ -2828,44 +2819,6 @@ public:
                            DenseMatrix &dshape) const;
 };
 
-/// Class for summation-by-parts operator on interval
-class SBP_SegmentElement : public NodalTensorFiniteElement
-{
-private:
-#ifndef MFEM_THREAD_SAFE
-   mutable Vector shape_x, dshape_x;
-#endif
-
-public:
-   SBP_SegmentElement(const int p);
-   virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
-   virtual void CalcDShape(const IntegrationPoint &ip,
-                           DenseMatrix &dshape) const;
-   
-   // ProjectDelta is leftover function from H1_SegmentElement
-   // virtual void ProjectDelta(int vertex, Vector &dofs) const;
-};
-
-/// Class for (diagonal-norm) summation-by-parts operator on triangles
-class SBP_TriangleElement : public NodalFiniteElement
-{
-private:
-#ifndef MFEM_THREAD_SAFE
-   mutable Vector shape_x, shape_y, shape_l, dshape_x, dshape_y, dshape_l, u;
-   mutable Vector ddshape_x, ddshape_y, ddshape_l;
-   mutable DenseMatrix du, ddu;
-#endif
-   DenseMatrix *Dx, *Dy;
-   DenseMatrixInverse Ti;
-
-public:
-   SBP_TriangleElement(const int p, const int Do);
-   virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
-   virtual void CalcDShape(const IntegrationPoint &ip,
-                           DenseMatrix &dshape) const;
-   virtual ~SBP_TriangleElement();
-};
-
 } // namespace mfem
 
 #endif

fem/fe_coll.cpp

@@ -274,10 +274,6 @@ FiniteElementCollection *FiniteElementCollection::New(const char *name)
          fec = new NURBSFECollection();
       }
    }
-   else if (!strncmp(name, "SBP_", 4))
-   {
-      fec = new C_SBPCollection(atoi(name+8), atoi(name+4));
-   }
    else
    {
       MFEM_ABORT("unknown FiniteElementCollection: " << name);
@@ -2546,159 +2542,4 @@ FiniteElementCollection *NURBSFECollection::GetTraceCollection() const
    return NULL;
 }
 
-C_SBPCollection::C_SBPCollection(const int p, const int dim)
-{
-   MFEM_VERIFY(p >= 0 && p <= 4, "C_SBPCollection requires 0 <= order <= 4.");
-   MFEM_VERIFY(dim == 2, "C_SBPCollection requires dim == 2.");
-
-   snprintf(c_SBPname, 32, "SBP_%dD_P%d", dim, p);
-
-   for (int g = 0; g < Geometry::NumGeom; g++)
-   {
-      C_SBPdof[g] = 0;
-      C_SBPElements[g] = NULL;
-   }
-   for (int i = 0; i < 2; i++)
-   {
-      SegDofOrd[i] = NULL;
-   }
-
-   C_SBPdof[Geometry::POINT] = 1;
-   C_SBPElements[Geometry::POINT] = new PointFiniteElement;
-
-   if (dim >= 1)
-   {
-      C_SBPdof[Geometry::SEGMENT] = p;
-
-      C_SBPElements[Geometry::SEGMENT] = new SBP_SegmentElement(p);
-
-      int nodeOrder0[] = {};
-      int nodeOrder1[1] = {0};
-      int nodeOrder2[2] = {0, 1};
-      int nodeOrder3[3] = {0, 1, 2};
-      int nodeOrder4[4] = {0, 1, 2, 3};
-
-      int revNodeOrder0[] = {};
-      int revNodeOrder1[1] = {0};
-      int revNodeOrder2[2] = {1, 0};
-      int revNodeOrder3[3] = {0, 2, 1};
-      int revNodeOrder4[4] = {1, 0, 3, 2};
-
-      switch (p)
-      {
-         case 0:
-            SegDofOrd[0] = new int[p];
-            SegDofOrd[1] = new int[p];
-            for (int i = 0; i < p; i++)
-            {
-               SegDofOrd[0][i] = nodeOrder0[i];
-               SegDofOrd[1][i] = revNodeOrder0[i];
-            }
-            break;
-         case 1:
-            SegDofOrd[0] = new int[p];
-            SegDofOrd[1] = new int[p];
-            for (int i = 0; i < p; i++)
-            {
-               SegDofOrd[0][i] = nodeOrder1[i];
-               SegDofOrd[1][i] = revNodeOrder1[i];
-            }
-            break;
-         case 2:
-            SegDofOrd[0] = new int[p];
-            SegDofOrd[1] = new int[p];
-            for (int i = 0; i < p; i++)
-            {
-               SegDofOrd[0][i] = nodeOrder2[i];
-               SegDofOrd[1][i] = revNodeOrder2[i];
-            }
-            break;
-         case 3:
-            SegDofOrd[0] = new int[p];
-            SegDofOrd[1] = new int[p];
-            for (int i = 0; i < p; i++)
-            {
-               SegDofOrd[0][i] = nodeOrder3[i];
-               SegDofOrd[1][i] = revNodeOrder3[i];
-            }
-            break;
-         case 4:
-            SegDofOrd[0] = new int[p];
-            SegDofOrd[1] = new int[p];
-            for (int i = 0; i < p; i++)
-            {
-               SegDofOrd[0][i] = nodeOrder4[i];
-               SegDofOrd[1][i] = revNodeOrder4[i];
-            }
-            break;
-         default:
-            mfem_error("SBP elements are currently only supported for 0 <= order <= 4");
-            break;
-
-      }
-   }
-
-   if (dim >= 2)
-   {
-      switch (p)
-      {
-         case 0:
-            C_SBPdof[Geometry::TRIANGLE] = 3 - 3 - 3*p;
-            break;
-         case 1:
-            C_SBPdof[Geometry::TRIANGLE] = 7 - 3 - 3*p;
-            break;
-         case 2:
-            C_SBPdof[Geometry::TRIANGLE] = 12 - 3 - 3*p;
-            break;
-         case 3:
-            C_SBPdof[Geometry::TRIANGLE] = 18 - 3 - 3*p;
-            break;
-         case 4:
-            C_SBPdof[Geometry::TRIANGLE] = 27 - 3 - 3*p;
-            break;
-         default:
-            mfem_error("SBP elements are currently only supported for 0 <= order <= 4");
-            break;
-      }
-
-      const int &TriDof = C_SBPdof[Geometry::TRIANGLE] + 3*C_SBPdof[Geometry::POINT] + 3*C_SBPdof[Geometry::SEGMENT];
-
-      C_SBPElements[Geometry::TRIANGLE] = new SBP_TriangleElement(p, TriDof);
-   }
-}
-
-const FiniteElement *C_SBPCollection::FiniteElementForGeometry(
-      Geometry::Type GeomType) const
-{
-   if (GeomType == Geometry::TRIANGLE || GeomType == Geometry::SEGMENT || GeomType == Geometry::POINT)
-   {
-
-   }
-   else
-   {
-      MFEM_ABORT("Unsupported geometry type " << GeomType);
-   }
-   return C_SBPElements[GeomType]; 
-}
-
-const int *C_SBPCollection::DofOrderForOrientation(Geometry::Type GeomType,
-                                                   int Or) const
-{
-   if (GeomType == Geometry::SEGMENT)
-   {
-      return (Or > 0) ? SegDofOrd[0] : SegDofOrd[1];
-   }
-   return NULL;
-}
-
-C_SBPCollection::~C_SBPCollection()
-{
-   delete [] SegDofOrd[0];
-   for (int g = 0; g < Geometry::NumGeom; g++)
-   {
-      delete C_SBPElements[g];
-   }
-}
-
 }

fem/fe_coll.hpp

@@ -884,30 +884,6 @@ public:
    virtual ~Local_FECollection() { delete Local_Element; }
 };
 
-/// Arbitrary order H1-conforming (continuous) Summation By Parts
-/// opperators.
-class C_SBPCollection : public FiniteElementCollection
-{
-
-protected:
-   char c_SBPname[32];
-   FiniteElement *C_SBPElements[Geometry::NumGeom];
-   int  C_SBPdof[Geometry::NumGeom];
-   int *SegDofOrd[2];
-
-public:
-   explicit C_SBPCollection(const int p, const int dim = 2);
-
-   virtual const FiniteElement *FiniteElementForGeometry(
-      Geometry::Type GeomType) const;
-   virtual int DofForGeometry(Geometry::Type GeomType) const
-   { return C_SBPdof[GeomType]; }
-   virtual const int *DofOrderForOrientation(Geometry::Type GeomType,
-                                             int Or) const;
-   virtual const char *Name() const { return c_SBPname; }
-   virtual ~C_SBPCollection();
-};
-
 }
 
 #endif

fem/gridfunc.cpp

@@ -2310,8 +2310,7 @@ double GridFunction::ComputeLpError(const double p, Coefficient &exsol,
       else
       {
          int intorder = 2*fe->GetOrder() + 1; // <----------
-         // ir = &(IntRules.Get(fe->GetGeomType(), intorder));
-         ir = &(IntRules.Get(*fe, intorder));
+         ir = &(IntRules.Get(fe->GetGeomType(), intorder));
       }
       GetValues(i, *ir, vals);
       T = fes->GetElementTransformation(i);

fem/intrules.cpp

@@ -907,30 +907,6 @@ const IntegrationRule &IntegrationRules::Get(int GeomType, int Order)
    return *(*ir_array)[Order];
 }
 
-/// Overloaded integration rule where FiniteElement is passed in, and the 
-/// OperatorType is used to determine which integration rules to use
-const IntegrationRule &IntegrationRules::Get(const FiniteElement &el, int Order)
-{
-   const IntegrationRule *ir = NULL;
-   switch (el.OperatorType)
-   {
-      case 0: // FE
-         ir = &Get(el.GetGeomType(), Order);
-         break;
-      case 1: // SBP
-         ir = &el.GetNodes(); // SBP type elements have collocated quadrature
-                               // notes and DOFs, weights are included in 
-                               // element construction so complete integration
-                               // rule is defined by the element's `Nodes`
-
-         break;
-      default:
-         MFEM_ABORT("Invalid OperatorType = " << el.OperatorType);
-         break;
-   }
-   return *ir;
-}
-
 void IntegrationRules::Set(int GeomType, int Order, IntegrationRule &IntRule)
 {
    Array<IntegrationRule *> *ir_array;

fem/intrules.hpp

@@ -18,20 +18,16 @@
 namespace mfem
 {
 
-class FiniteElement; // Forward declaration of class to alow overloaded Get()
-
 /* Classes for IntegrationPoint, IntegrationRule, and container class
    IntegrationRules.  Declares the global variable IntRules */
 
 /// Class for integration point with weight
 class IntegrationPoint
 {
-private:
-   int index;
 public:
    double x, y, z, weight;
 
-   void Init(int idx) { x = y = z = weight = 0.0; index = idx;}
+   void Init() { x = y = z = weight = 0.0; }
 
    void Set(const double *p, const int dim)
    {
@@ -83,8 +79,6 @@ public:
    void Set1w(const double x1, const double w) { x = x1; weight = w; }
 
    void Set1w(const double *p) { x = p[0]; weight = p[1]; }
-
-   const int GetIdx() const {return index;}
 };
 
 /// Class for an integration rule - an Array of IntegrationPoint.
@@ -218,7 +212,7 @@ public:
    {
       for (int i = 0; i < this->Size(); i++)
       {
-         (*this)[i].Init(i);
+         (*this)[i].Init();
       }
    }
 
@@ -354,9 +348,6 @@ public:
    /// Returns an integration rule for given GeomType and Order.
    const IntegrationRule &Get(int GeomType, int Order);
 
-   /// Returns an integration rule for an element and order
-   const IntegrationRule &Get(const FiniteElement &el, int Order);
-
    void Set(int GeomType, int Order, IntegrationRule &IntRule);
 
    void SetOwnRules(int o) { own_rules = o; }

fem/lininteg.cpp

@@ -38,8 +38,7 @@ void DomainLFIntegrator::AssembleRHSElementVect(const FiniteElement &el,
    {
       // ir = &IntRules.Get(el.GetGeomType(),
       //                    oa * el.GetOrder() + ob + Tr.OrderW());
-      // ir = &IntRules.Get(el.GetGeomType(), oa * el.GetOrder() + ob);
-      ir = &IntRules.Get(el, oa * el.GetOrder() + ob); // New overload
+      ir = &IntRules.Get(el.GetGeomType(), oa * el.GetOrder() + ob);
    }
 
    for (int i = 0; i < ir->GetNPoints(); i++)

linalg/hypre.hpp

@@ -468,16 +468,13 @@ public:
    { internal::hypre_ParCSRMatrixSetConstantValues(A, value); return *this; }
 
    /** Perform the operation `*this += B`, assuming that both matrices use the
-       same row and column partitions and the same col_map_offd arrays, or B has
-       an empty off-diagonal block. We also assume that the sparsity pattern of
-       `*this` contains that of `B`. */
+       same row and column partitions and the same col_map_offd arrays. We also
+       assume that the sparsity pattern of `*this` contains that of `B`. */
    HypreParMatrix &operator+=(const HypreParMatrix &B) { return Add(1.0, B); }
 
    /** Perform the operation `*this += beta*B`, assuming that both matrices use
-       the same row and column partitions and the same col_map_offd arrays, or
-       B has an empty off-diagonal block. We also assume that the sparsity
-       pattern of `*this` contains that of `B`. For a more general case consider
-       the stand-alone function ParAdd described below. */
+       the same row and column partitions and the same col_map_offd arrays. We
+       also assume that the sparsity pattern of `*this` contains that of `B`. */
    HypreParMatrix &Add(const double beta, const HypreParMatrix &B)
    {
       MFEM_VERIFY(internal::hypre_ParCSRMatrixSum(A, beta, B.A) == 0,

linalg/hypre_parcsr.hpp

@@ -115,9 +115,8 @@ hypre_ParCSRMatrixAdd(hypre_ParCSRMatrix *A,
                       hypre_ParCSRMatrix *B);
 
 /** Perform the operation A += beta*B, assuming that both matrices use the same
-    row and column partitions and the same col_map_offd arrays, or B has an empty
-    off-diagonal block. We also assume that the sparsity pattern of A contains
-    that of B. */
+    row and column partitions and the same col_map_offd arrays. We also assume
+    that the sparsity pattern of A contains that of B. */
 HYPRE_Int
 hypre_ParCSRMatrixSum(hypre_ParCSRMatrix *A,
                       HYPRE_Complex       beta,

linalg/ode.cpp

@@ -561,6 +561,102 @@ void GeneralizedAlphaSolver::Step(Vector &x, double &t, double &dt)
 }
 
 
+void IMEX_BE_FE::Init(TimeDependentOperator &_f)
+{
+   ODESolver::Init(_f);
+   k_imp.SetSize(f->Width());
+   y.SetSize(f->Width());
+   k_exp.SetSize(f->Width());
+}
+
+void IMEX_BE_FE::Step(Vector &x, double &t, double &dt)
+{
+   f->ExplicitMult(x, k_exp);
+   add(x, dt, k_exp, y);
+
+   f->SetTime(t + dt);
+   f->ImplicitSolve(dt, y, k_imp);
+
+   x.Add(dt, k_exp);
+   x.Add(dt, k_imp);
+   t += dt;
+}
+
+
+void IMEXRK2::Init(TimeDependentOperator &_f)
+{
+   ODESolver::Init(_f);
+   f = ODESolver::f;
+   k_imp.SetSize(f->Width());
+   k_exp.SetSize(f->Width());
+   y.SetSize(f->Width());
+   z.SetSize(f->Width());
+}
+
+void IMEXRK2::Step(Vector &x, double &t, double &dt)
+{
+   double gamma = 1 - sqrt(2)/2;
+   double delta = -2*sqrt(2)/3;
+
+   // The method is given by
+   // k1_exp = f(u)
+   // k1_imp = g(u + gamma*dt*k1_exp + gamma*dt*k1_imp)
+   // k2_exp = f(u + gamma*dt*k1_exp + gamma*dt*k1_imp)
+   // k2_imp = g(u + delta*dt*k1_exp + (1-gamma)*dt*k1_imp
+   //              + (1-delta)*dt*k2_exp + gamma*dt*k2_imp)
+   // k3_exp = f(u + delta*dt*k1_exp + (1-gamma)*dt*k1_imp
+   //              + (1-delta)*dt*k2_exp + gamma*dt*k2_imp)
+   // u_new = u + dt*((1-gamma)*k1_imp + (1-gamma)*k2_exp
+   //                 + gamma*k2_imp + gamma*k3_exp)
+
+   // Take first explicit step
+   // k1_exp = f(u)
+   f->ExplicitMult(x, k_exp);
+   // b corresponding to this stage is zero, so don't add to solution
+
+   // Solve first implicit step
+   // y = u + gamma*dt*k1_exp
+   add(x, gamma*dt, k_exp, y);
+   // Solve x1_imp = g(u + gamma*dt*k1_exp + gamma*dt*k1_imp)
+   f->SetTime(t + gamma*dt);
+   f->ImplicitSolve(gamma*dt, y, k_imp);
+   // x = u + (1-gamma)*dt*k1_imp
+   x.Add((1-gamma)*dt, k_imp);
+
+   // Begin setting up rhs for second solve
+   // z = u + (1-gamma)*dt*k_imp + delta*dt*k_exp
+   add(x, delta*dt, k_exp, z);
+
+   // Take second explicit step
+   // y = x + gamma*dt*k1_exp + gamma*dt*k1_imp
+   y.Add(gamma*dt, k_imp);
+   // k2_exp = f(x + gamma*dt*k1_exp + gamma*dt*k1_imp)
+   f->ExplicitMult(y, k_exp);
+   // x = u + (1-gamma)*dt*k1_imp + (1-gamma)*dt*k2_exp
+   x.Add((1-gamma)*dt, k_exp);
+
+   // Finish formoing rhs
+   // z = x + (1-gamma)*dt*k1_imp + delta*dt*k1_exp + (1-delta)*dt*k2_exp
+   z.Add((1-delta)*dt, k_exp);
+
+   // Solve second implicit step for k2_imp
+   f->SetTime(t + dt);
+   f->ImplicitSolve(gamma*dt, z, k_imp);
+   // x = u + (1-gamma)*dt*k1_imp + (1-gamma)*dt*k2_exp + gamma*dt*k2_imp
+   x.Add(gamma*dt, k_imp);
+
+   // Take final explicit step for k3_exp
+   z.Add(gamma*dt, k_imp);
+   f->ExplicitMult(z, k_exp);
+
+   // x = u + (1-gamma)*dt*k1_imp + (1-gamma)*dt*k2_exp + gamma*dt*k2_imp
+   //       + gamma*dt*k3_exp
+   x.Add(gamma*dt, k_exp);
+
+   t += dt;
+}
+
+
 void
 SIASolver::Init(Operator &P, TimeDependentOperator & F)
 {

linalg/ode.hpp

@@ -305,6 +305,33 @@ public:
 };
 
 
+/// IMEX Backward-Forward Euler ODE solver
+class IMEX_BE_FE : public ODESolver
+{
+protected:
+   Vector k_exp, k_imp, y;
+
+public:
+   virtual void Init(TimeDependentOperator &_f);
+
+   virtual void Step(Vector &x, double &t, double &dt);
+};
+
+/** Second-order IMEX (2,3,2) method, from "Implicit-explicit Runge-Kutta
+    methods for time-dependent partial differential equations" by Ascher, Ruuth
+    and Spiteri, Applied Numerical Mathematics (1997). */
+class IMEXRK2 : public ODESolver
+{
+protected:
+   Vector k_exp, k_imp, y, z;
+
+public:
+   virtual void Init(TimeDependentOperator &_f);
+
+   virtual void Step(Vector &x, double &t, double &dt);
+};
+
+
 /// The SIASolver class is based on the Symplectic Integration Algorithm
 /// described in "A Symplectic Integration Algorithm for Separable Hamiltonian
 /// Functions" by J. Candy and W. Rozmus, Journal of Computational Physics,