Cleaned up unordered_map implementation and added safer float comparison as a fallback for out_of_range errors thrown by the map.

Merging sbp-dev into hashtable-sbp-dev to test hashtable implementation with more recent sbp-dev branch before merging back to sbp-dev

.gitignore

@@ -39,6 +39,9 @@ doc/CodeDocumentation
 *.dSYM
 .DS_Store
 
+# Editor files
+.vscode
+
 # Example and miniapp binaries and outputs
 
 examples/ex[1-9]
@@ -85,6 +88,7 @@ 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,7 +28,6 @@ list(APPEND ALL_EXE_SRCS
   ex19.cpp
   ex20.cpp
   ex22.cpp
-  ex23.cpp
   )
 
 if (MFEM_USE_MPI)
@@ -54,7 +53,6 @@ if (MFEM_USE_MPI)
     ex19p.cpp
     ex20p.cpp
     ex22p.cpp
-    ex23p.cpp
     )
 endif()
 

examples/ex23.cpp

@@ -1,734 +0,0 @@
-//                               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

@@ -1,797 +0,0 @@
-//                       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

@@ -0,0 +1,359 @@
+//                                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>
+#include <chrono>
+
+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 = 1;
+	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, 1000, 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 ex23
+  ex18 ex19 ex20 ex22 exSBP
 PAR_EXAMPLES = ex1p ex2p ex3p ex4p ex5p ex6p ex7p ex8p ex9p ex10p ex11p ex12p\
- ex13p ex14p ex15p ex16p ex17p ex18p ex19p ex20p ex22p ex23p
+ ex13p ex14p ex15p ex16p ex17p ex18p ex19p ex20p ex22p
 
 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* Example23*
+	@rm -rf Example5* Example9* Example15* Example16*
 	@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,4 +126,3 @@ 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

@@ -415,6 +415,13 @@ void DiffusionIntegrator::AssembleElementMatrix
       {
          ir = &RefinedIntRules.Get(el.GetGeomType(), order);
       }
+      else if (el.Space() == FunctionSpace::SBPk)
+      {
+         ir = &el.GetNodes(); // SBP elements have collocated quadrature nodes
+                              // and DOFs, weights are included in element
+                              // construction so complete integration rule is
+                              // defined by the element's `Nodes`
+      }
       else
       {
          ir = &IntRules.Get(el.GetGeomType(), order);

fem/fe.cpp

@@ -11845,6 +11845,376 @@ 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::SBPk)
+{
+   /// Header file including SBP Dx and Dy matrix data
+   #include "fe_sbp.hpp"
+
+   // Create Dx and Dy matrixes
+   Dx = new DenseMatrix(Dof);
+   Dy = new DenseMatrix(Dof);
+   
+   // 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;
+   }
+
+   // populate unordered_map with mapping from IntPoint address to index
+   for (int i = 0; i < Dof; i++)
+   {
+      ipIdxMap[&(Nodes.IntPoint(i))] = i;
+   }
+}
+
+/// 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
+{
+   int ipIdx;
+   try
+   {
+      ipIdx = ipIdxMap.at(&ip);
+   }
+   catch (const std::out_of_range& oor)
+   // error handling code to handle cases where the pointer to ip is not
+   // in the map. Problems arise in GridFunction::SaveVTK() ->  GridFunction::GetValues()
+   // which calls CalcShape() with an `IntegrationPoint` defined by a refined
+   // geometry type. Since the IntegrationPoint is not in Nodes, its address is
+   // not in the ipIdxMap, and an out_of_range error is thrown. This code catches
+   // the error and uses float comparisons to determine the IntegrationPoint
+   // index.
+   {
+      double tol = 1e-12;
+      for (int i = 0; i < Dof; i++)
+      {
+         double delta_x = ip.x - Nodes.IntPoint(i).x;
+         double delta_y = ip.y - Nodes.IntPoint(i).y;
+         if (delta_x*delta_x + delta_y*delta_y < tol)
+         {
+            ipIdx = i;
+            break;
+         }
+      }
+   }
+   shape = 0.0;
+   shape(ipIdx) = 1.0;
+}
+
+/// 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;
+   try
+   {
+      ipIdx = ipIdxMap.at(&ip);
+   }
+   catch (const std::out_of_range& oor)
+   // error handling code to handle cases where the pointer to ip is not
+   // in the map. Problems arise in GridFunction::SaveVTK() ->  GridFunction::GetValues()
+   // which calls CalcShape() with an `IntegrationPoint` defined by a refined
+   // geometry type. Since the IntegrationPoint is not in Nodes, its address is
+   // not in the ipIdxMap, and an out_of_range error is thrown. This code catches 
+   // the error and uses float comparisons to determine the IntegrationPoint
+   // index.
+   {
+      double tol = 1e-12;
+      for (int i = 0; i < Dof; i++)
+      {
+         double delta_x = ip.x - Nodes.IntPoint(i).x;
+         double delta_y = ip.y - Nodes.IntPoint(i).y;
+         if (delta_x*delta_x + delta_y*delta_y < tol)
+         {
+            ipIdx = i;
+            break;
+         }
+      }
+   }
+   dshape = 0.0;
+
+   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()
+   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

@@ -19,6 +19,7 @@
 #include "geom.hpp"
 
 #include <map>
+#include <unordered_map>
 
 namespace mfem
 {
@@ -126,7 +127,8 @@ public:
    {
       Pk, ///< Polynomials of order k
       Qk, ///< Tensor products of polynomials of order k
-      rQk ///< Refined tensor products of polynomials of order k
+      rQk,///< Refined tensor products of polynomials of order k
+      SBPk///< Summation-by-parts operator of order k with no explicit basis
    };
 };
 
@@ -2819,6 +2821,44 @@ 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;
+   std::unordered_map<const IntegrationPoint*, int> ipIdxMap;
+
+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,6 +274,10 @@ 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);
@@ -2542,4 +2546,159 @@ 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,6 +884,30 @@ 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/fe_sbp.hpp

@@ -0,0 +1,141 @@
+/// SBP Triangle Dx and Dy matrix data arrays, stored column major format
+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};

fem/gridfunc.cpp

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

fem/lininteg.cpp

@@ -36,9 +36,16 @@ void DomainLFIntegrator::AssembleRHSElementVect(const FiniteElement &el,
    const IntegrationRule *ir = IntRule;
    if (ir == NULL)
    {
-      // ir = &IntRules.Get(el.GetGeomType(),
-      //                    oa * el.GetOrder() + ob + Tr.OrderW());
-      ir = &IntRules.Get(el.GetGeomType(), oa * el.GetOrder() + ob);
+      if (el.Space() == FunctionSpace::SBPk)
+      {
+         ir = &el.GetNodes();
+      }
+      else 
+      {
+         // ir = &IntRules.Get(el.GetGeomType(),
+         //                    oa * el.GetOrder() + ob + Tr.OrderW());
+         ir = &IntRules.Get(el.GetGeomType(), oa * el.GetOrder() + ob);
+      }
    }
 
    for (int i = 0; i < ir->GetNPoints(); i++)

linalg/hypre.hpp

@@ -468,13 +468,16 @@ 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. 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, or B has
+       an empty off-diagonal block. 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. We
-       also assume that the sparsity pattern of `*this` contains that of `B`. */
+       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. */
    HypreParMatrix &Add(const double beta, const HypreParMatrix &B)
    {
       MFEM_VERIFY(internal::hypre_ParCSRMatrixSum(A, beta, B.A) == 0,

linalg/hypre_parcsr.hpp

@@ -115,8 +115,9 @@ 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. We also assume
-    that the sparsity pattern of A contains that of B. */
+    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. */
 HYPRE_Int
 hypre_ParCSRMatrixSum(hypre_ParCSRMatrix *A,
                       HYPRE_Complex       beta,

linalg/ode.cpp

@@ -561,102 +561,6 @@ 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,33 +305,6 @@ 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,